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Abstract 


A  numerical  method  has  been  developed  for  solving  the  complete  compressible  Navier-Stokes 
equations.  The  method  is  applicable  for  Direct  Numerical  Simulations  (DNS)  and  Large-Eddy 
Simulations  (LES)  and  was  used  here  to  study  the  evolution  of  three-dimensional  disturbances 
in  the  laminar  and  turbulent  near  wake  of  axisymmetric  bluff  bodies  with  a  blunt  base  in 
supersonic  flows.  The  main  objective  of  this  research  is  to  investigate  the  time  dependent 
behavior  of  these  distmrbances  and  their  influence  on  and  interaction  with  the  global  flow  field. 
The  equations  are  solved  in  a  cylindrical  coordinate  system  using  finite  difference  approximations 
of  fourth-order  accuracy  in  axial  and  radial  directions  and  and  a  fourth-order  accurate  explicit 
Runge-Kutta  scheme  for  the  time  integration.  A  pseudo-spectral  method  is  employed  in  the 
azimuthal  direction.  Direct  Numerical  Simulations  (DNS)  were  performed  for  a  subsonic  free 
stream  Mach  number  of  Moo  =  0.2  and  for  supersonic  free  stream  Mach  numbers  of  Moo  =  1-2 
and  Moo  =  2.46.  Large-Eddy  Simulations  (LES)  were  carried  out  for  a  subsonic  free  stream 
Mach  number  of  Moo  =  0.2  and  a  global  Reynolds  number  oi  Brd  =  2,000  and  for  a  supersonic 
free  stream  Mach  number  of  Moo  =  2.46  and  global  Reynolds  numbers  of  Bed  =  30,000  and 
Bed  =  100,000.  Comparison  of  the  instantaneous  flow  field  for  subsonic  calculations  with 
water  channel  experiments  and  incompressible  simulations  show  good  qualitative  agreement. 
An  absolute  instability  with  regard  to  helical  disturbances  was  found  for  the  subsonic  flow  at 
Bed  =  1,000  and  for  the  supersonic  flows  for  =  1.2  and  Bed  >  4,000  and  for  Mqo  =  2.46 
and  Bed  >  30, 000.  Small  disturbances  appear  in  the  flow  field  near  the  corner  of  the  base.  As 
the  disturbances  are  propagating  downstream  they  grow  and  form  intense  vortical  structures. 
These  structures  have  a  strong  influence  on  the  flow  field,  which  results  in  a  drastic  change  of 
the  base  pressure  distribution  and  thus  of  the  base  drag. 
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NOMENCLATURE 

Roman  Symbols 
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speed  of  sound 

Cp 

specific  heat  at  constant  pressme 

Cs 

Smagorinsky  constant  for  Large-Eddy  Simulations 
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specific  heat  at  constant  volume 
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diameter  of  the  axisymmetric  body 
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total  energy  per  unit  mass 

Er 

resolved  scale  total  energy 
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frequency 

k 

Fourier  mode  of  Fourier  decomposition  in  azimuthal  direction 
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0'  unfiltered  part  for  spatial  filter 

0"  unfiltered  part  for  Favre  filter 

Subscripts 

oo  dimensional  quantities  evaluated  at  the  free  stream 

X,  r,  6  subscripts  indicating  the  component  of  a  vector  or  tensor 

Other  Symbols 

0  filtered  quantities 

0  Favre  filtered  (mass-averaged)  quantities 

<  >  root-mean-squared  quantity  (RMS) 


Chapter  1 


Introduction 


The  aerodynamic  drag  of  bodies  of  revolution  with  a  blunt  base  is  greatly  affected  by  the 
low  pressure  immediately  downstream  of  the  base.  Especially  for  the  case  of  supersonic  flight, 
even  small  changes  in  the  flow  behavior  of  the  wake  may  affect  the  performance  of  the  entire 
flight  vehicle,  e.g.,  missiles,  rockets,  or  projectiles.  Flight  tests  with  projectiles  (U.S.  Army  M549 
projectile)  have  shown  that  for  supersonic  free  stream  Mach  numbers  the  fraction  of  aerodynamic 
drag  due  to  the  low  base  pressure  can  be  as  high  as  35  percent  (compare  flgure  1.1)  of  the  total 
drag.  This  suggests  that  controlling  the  near-wake  flow  can  have  a  signiflcant  influence  on  the 
aerodynamic  drag  of  bluff  bodies  with  a  blunt  base. 


Macfa  Number 

Figure  1.1:  Total  drag  and  base  drag  for  axisymmetric  bluff  bodies  with  a  blunt  base  in  free 
flight  experiments  [from  Rollstin  (1987)]. 

For  this  reason,  numerous  research  efforts,  experimental,  theoretical,  and  numerical,  have 
focused  on  investigating  the  flow  fleld  in  the  near  wake  of  blunt  bodies.  At  the  early  stages 
the  main  issue  of  these  investigations  has  been  the  prediction  of  the  base  drag  for  different 
geometries  in  order  to  provide  necessary  information  for  the  design  of  new  aerodynamic  objects 
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(such  as  projectiles,  rockets,  supersonic  planes,  etc.).  In  more  recent  years,  several  techniques  for 
reducing  the  base  drag  have  been  suggested  and  investigated,  such  as  boat-tailing,  base  bleed, 
base  biurning,  base  combustion,  etc.  [see,  for  example,  SaJiu  et  al.(1985),  Nietubicz  and  Gibeling 
(1993),  Sahu  and  Heavey  (1995),  and  Mathur  and  Dutton  (1996a, b)].  It  has  been  found,  that 
with  all  these  techniques  the  base  drag  can  be  reduced  by  a  considerable  amount. 


Figure  1.2:  Schematic  of  the  flow  fleld  of  the  neair  wcike  of  axisymmetric  blufll  bodies  with  a 
blunt  base  aligned  with  a  supersonic  free  stream  [from  Herrin  &  Dutton  (1994c)]. 

A  schematic  diagram  of  the  time-averaged  supersonic  flow  field  of  the  near  wake  of  an 
axisymmetric  blufi'body  with  a  blunt  base  is  shown  in  figure  1.2.  A  boundary  layer  separates  at 
the  bcise  corner  and  forms  an  axisymmetric  shear  layer.  Outside  the  subsonic  layer,  expansion 
waves  occur  that  cause  the  free  stream  to  bend  towards  the  axis.  Further  downstream  the 
free  stream  is  realigned  by  recompression  waves.  A  region  of  recirculation,  associated  with 
low  pressure,  is  formed  immediately  downstream  of  the  base.  Typically  the  length  of  this 
recirculation  region  is  on  the  order  of  one  to  three  base  diameters,  depending  on  the  global 
Reynolds  number  (ife/j)  and  the  free  stream  Mach  number  {M^). 

1.1  Technical  Background 

1.1.1  Investigations  of  the  Time- Averaged  Flow  Field 

At  first,  for  investigating  the  fundamental  physics  of  the  near- wake  behavior  experiments  were 
used  almost  exclusively.  However,  wind  tunnel  experiments  suffer  from  the  diflficulty  of  properly 
supporting  the  (axisymmetric)  base  model  so  as  not  to  cause  undue  effects  on  the  flow  field 
behind  the  base.  [Chapman  (1951)]  conducted  experiments  mecisuring  the  base  pressure  behind 
a  backward  facing  step  and  a  circular  cylinder  for  supersonic  flows.  He  found  that  the  support 
of  the  model,  which  in  his  case  was  supported  from  the  rear,  in  general  has  a  large  efiect  on  the 
base  pressure  if  the  sting  diameter  is  larger  than  30  percent  of  the  model  diameter.  Another 
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experimental  study  of  the  effect  of  rear  sting  supports  on  the  base  pressure  was  conducted 
by  [Donaldson  (1955)].  He  observed  that  the  base  pressure  dropped  monotonically  with  an 
increasing  sting  diameter  up  to  80  percent  of  the  base  diameter.  Even  a  support  rod  with  a 
diameter  of  20  percent  of  the  base  diameter  had  a  significant  effect  on  the  base  pressure,  in 
contrast  to  the  findings  of  [Chapman  (1951)]. 

Because  of  the  problems  arrising  from  the  rear  support  of  the  model,  several  other  methods 
of  support  have  been  investigated  in  the  following  years.  For  example,  [Dayman  (1963)]  investi¬ 
gated  the  effect  of  supporting  wires  on  the  wakes  of  spheres  and  cones.  He  concluded  that  any 
wire  capable  of  supporting  the  model  has  a  significant  effect  on  the  flow  field.  Thus,  the  only 
support  that  might  not  significantly  influence  the  wake  flow  field  is  from  the  upstream  position 
of  the  body,  which  has  been  used  in  most  of  the  experimental  research  efforts  in  the  following 
years  [see,  for  example,  [Dutton  &  Addy  (1993)],  Herrin  &  Dutton  (1991,  1994a-c),  and  Mathur 
&  Dutton  (1995,  1996a,b)].  In  addition,  local  pressure  and  velocity  measurements  in  supersonic 
wakes  have  proven  to  be  very  difficult.  The  pressiure  along  the  body  can  be  measured  by  pressure 
tabs,  but  measurements  within  the  flow  field  are  very  difficult,  because  a  probe  has  to  be  inserted 
into  the  flow  and  thus  the  flow  field  is  altered.  That  is  the  reason  why  different,  nonintruding 
methods  such  as  laser  Doppler  velocimetry  have  been  developed.  In  many  recent  experimen¬ 
tal  studies  [see,  for  example,  [Avidor  &  Schneiderman  (1975)],  [Herrin  &  Dutton  (1991)]]  near 
axisymmetric  wakes  have  been  investigated  using  this  method. 

In  recent  years,  many  experimental  investigations  also  focused  on  modifying  the  base  flow, 
in  particular  the  base  pressure,  using  several  different  techniques  such  as  boat-tailing,  base 
bleed,  base  burning,  and  base  combustion.  Valentine  &  Przirembel  (1970),  for  example,  in¬ 
vestigated  the  effect  of  base  injection  on  a  turbulent  axisymmetric  wake  at  Mach  4.  The  ef¬ 
fect  of  base  injection  on  the  base  pressure  was  dependent  on  the  injection  technique.  They 
observed  an  increase  in  the  base  pressure  ratio  of  nearly  100  percent.  Another  way  of  in¬ 
fluencing  the  base  pressure  was  assessed  by  [Neale  et  al.(1978)].  They  studied  the  effect  of 
axisymmetric  external  compression  fields  simulating  external  burning  on  the  base  pressure  for 
a  cylindrical  body  in  a  supersonic  flow  at  Mach  3.  Their  conclusion  was  that  the  base  drag 
could  be  completely  eliminated  by  external  burning.  There  have  been  many  more  research 
efforts  investigating  the  effects  of  these  techniques  to  modify  the  base  flow  [see,  for  exam¬ 
ple,  [Cortwright  &  Schroeder  (1951)],  [Reid  &  Hastings  (1959)],  [Clayden  &  Bowman  (1968)], 
[Bowman  &  Clayden  (1968)],  [Hubbartt  et  al.  (1981)],  Ding  et  al.  (1992)].  These  experiments 
have  shown  that  there  are  potentially  considerable  rewards  for  altering  the  flow  field.  However, 
it  is  still  not  understood  why  certain  measures  are  more  effective  than  others  and  what  the 
optimal  parameters  should  be.  The  reason  for  this  is  that  the  fundamental  physical  phenomena 
that  are  responsible  for  the  effect  of  the  base  flow  modifications  are  not  yet  understood. 

More  recently,  Dutton  and  his  coworkers  started  comprehensive  experimental  investigation 
of  the  near  wake  of  an  axisymmetric  body  with  a  blunt  base  in  a  supersonic  flow  at  Mach  2.46 
[see,  for  example,  Herrin  &:  Dutton  (1991,  1994a-c),  and  Mathur  &  Dutton  (1995,  1996a,b)]. 
In  their  research,  they  have  used  pressure  taps  along  the  base  to  measure  the  base  pressure 
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and  spark-schlieren  photography  and  oil-streak  visualization  to  qualitatively  assess  the  flow 
field.  Special  emphasis  was  placed  on  designing  the  experimental  facilities  such  that  a  highly 
uniform,  axisymmetric  flow  field  could  be  obtained  without  any  interference  of  the  support 
and  the  wind  tunnel  walls  on  the  wake  region.  It  was  found  that  the  base  pressure  is  virtually 
constant  along  the  base.  In  addition  to  the  base  pressure  measurements  they  used  Laser  Doppler 
Velocimetry  (LDV)  to  determine  the  time-averaged  flow  field  and  rms  values  of  the  turbulent 
fluctuations.  The  measmrements  of  turbulent  fluctuations  in  the  flow  field  revealed  that  the 
highest  turbulence  levels  occur  within  the  free  shear  layer  that  separates  the  recirculating  region 
from  the  free  stream.  Dutton  and  his  co-workers  have  also  investigated  the  efiFects  of  different 
techniques  of  drag  reduction,  including  base  bleed  and  boat-tailing.  Their  results  show  that 
all  these  techniques  can  significantly  lower  the  base  drag.  At  the  same  time,  all  investigated 
measures  of  drag  reduction  also  lowered  the  level  of  the  turbulent  fluctuations  in  the  free  shear 
layer.  This  indicates  that  there  might  be  a  connection  between  the  turbulent  energy  within  the 
shear  layer  and  the  pressmre  drag  of  the  aerodynamic  body. 

In  the  same  publication  where  [Chapman  (1951)]  discussed  his  experimental  measurements 
of  compressible  wakes,  he  also  presented  a  semi-empirical  theory  for  calculating  the  flow  field  of 
laminar,  viscous  wake  flow.  Chapman  had  concluded  that  inviscid  theory  was  not  insufiicient  to 
accurately  model  this  flow  field.  He  was  able  to  predict  the  base  pressure  for  two  dimensional 
wakes  with  reasonable  accmacy,  however,  he  did  not  succeed  in  predicting  the  base  pressure 
for  axisymmetric  wakes  without  rear  sting  support.  Another  semi-empirical  theory  for  the 
prediction  of  turbulent  flows  over  a  two  dimensional  backward  facing  step  and  an  axisymmetric 
body  with  rear  sting  support  wcis  presented  by  [Korst  (1956)].  However,  the  theory  was  also  not 
applicable  for  bodies  without  a  rear  sting.  A  similar  theory  was  developed  by  [McDonald  (1965)] 
which  also  gave  fair  agreement  with  experimental  results  for  the  base  pressure  of  axisymmtric 
bodies  with  a  rear  support.  A  newer  theory,  developed  by  [Mueller  (1968)],  was  able  to  calculate 
the  base  pressure  for  axisymmetric  bodies  without  rear  support  fairly  well.  But  he  still  had  to 
rely  on  experimental  data  to  find  the  recompression  location  for  a  given  shape,  which  made  it 
impossible  to  predict  the  base  pressure  for  a  new  geometry.  Several  theories  have  been  developed 
later  trying  to  predict  the  recompression  location  (e.g.  [Addy  (1970)]),  but  they  still  had  to  rely 
on  experimental  data. 

With  the  development  of  increasingly  powerful  supercomputers,  numerical  simulations  of 
complex  flows,  such  as  the  axisymmetric  waJce,  have  become  more  and  more  feasible.  The  big 
advantage  of  numerical  simulations  is  that  problems  associated  with  wind  tunnel  interference, 
model  support,  probe  intrusion,  etc.,  are  not  present.  In  addition  to  possibly  providing  fur¬ 
ther  understanding  of  the  relevant  physics  involved,  these  calculations  were  motivated  by  the 
considerable  challenge  that  this  complicated  flow  field  provides  for  computational  fluid  dynami- 
cists.  The  computational  challenge  arises  mainly  from  the  combination  of  shock  waves,  thin  free 
shear  layers,  boundary  layers,  and  recirculating  regions;  associated  with  this  are  highly  disparate 
length  scales  and  local  regions  of  very  high  gradients.  Reliable  and  realistic  computations  can 
therefore  only  be  performed  if  the  high  gradients  can  be  adequately  resolved.  Furthermore, 


CHAPTER  1.  INTRODUCTION 


7 


supersonic  wake  flows  in  almost  all  practical  applications  are  turbulent,  requiring  adequate  tur¬ 
bulence  models.  In  practically  all  previous  numerical  efforts,  only  the  steady  flow  fleld  was 
calculated.  Those  computations  were  based  on  the  Reynolds-Averaged  Navier-Stokes  equations 
(RANS)  employing  typically  algebraic  or  one-  and  two-equation  turbulence  models,  such  as,  for 
example,  the  algebraic  eddy  viscosity  model  of  Baldwin  and  Lomax,  the  two-equation  k  —  e 
model,  or  the  k  —  u  model. 

Numerous  research  efforts  on  calculating  turbulent  supersonic  wake  flows  have  been  reported 
in  the  literature.  The  Advisory  Group  for  Aerospace  Research  and  Development  (AGARD) 
formed  a  group  to  review  methods  to  predict  nozzle  afterbody  flows  in  1982.  A  summary  of  the 
results  is  given  by  Putnam  &  Bissinger  (1985).  It  was  concluded  that  at  the  time  Navier-Stokes 
simulations  were  able  to  predict  flow  flelds  accurately  up  to  the  point  of  separation,  but  then 
became  unreliable.  Since  then  there  has  been  an  increased  interest  in  predicting  wake  flows  using 
numerical  techniques.  [Sullins  et  al.  (1982)],  for  example,  calculated  the  turbulent  mean  flow  of 
a  two  dimensional  compressible  wake  with  and  without  base  injection.  For  the  determination 
of  the  turbulent  flow  they  used  a  zero  equation  relaxation  eddy  viscosity  model.  Comparing 
the  pressure  along  the  centerline  of  the  wake  showed  fair  agreement  with  experimental  results. 
Unfortunately  no  comparison  for  the  pressure  distribution  over  the  base  was  provided. 

The  flow  of  an  axisymmetric  afterbody  at  a  free  stream  Mach  number  of  1.343  with  and 
without  a  centered  jet  at  the  base  was  investigated  by  Sahu  &:  Nietubicz  (1984)  using  a  thin 
layer  approximation  of  the  Navier-Stokes  equation.  They  predicted  streamlines  for  the  near  wake 
region  and  also  the  pressure  distribution  over  the  base.  They  found  that  the  centered  jet  reduces 
the  base  pressure  for  the  case  of  a  high  velocity  jet.  No  comparison  with  experimental  data  was 
provided.  A  detailed  veriflcation  of  Navier-Stokes  calculations  was  done  by  Nietubicz  &  Sturek 
(1988).  The  results  show  very  good  agreement  with  experiments  for  most  flow  parameters  in¬ 
cluding  the  base  drag.  The  study  showed  that  current  Navier-Stokes  codes  axe  able  to  calculate 
supersonic  wake  flows  and  find  agreement  with  experimental  observations  in  terms  of  certain 
global  features  fairly  well.  However,  in  all  cases  certain  parameters  for  the  employed  turbulence 
models  have  to  be  set  according  to  experimental  results.  Numerical  computations  of  supersonic 
base  flow  with  special  emphasis  on  tmrbulence  modeling  were  presented  by  [Sahu  (1992)].  Com¬ 
paring  the  results  obtained  by  a  thin  layer  approximation  of  the  Navier-Stokes  equations  using 
two  algebraic  eddy  viscosity  models  and  a  two-equation  k-e  model,  he  found  that  both  algebraic 
tmbulence  models  predicted  an  incorrect  base  pressure,  while  the  k-e  model  showed  very  good 
agreement  with  experimental  results.  However,  the  turbulence  energy  distribution  in  the  wake 
did  not  agree  with  the  experimental  results  of  [Herrin  &  Dutton  (1991)]. 

In  spite  of  the  these  shortcomings  of  the  RANS  calculations,  recent  applications  have  at¬ 
tacked  increasingly  difficult  situations.  As  control  of  wake  flows  is  now  being  considered  as  a 
means  of  drag  reduction,  numerical  simulations  have  been  performed  to  investigate,  for  exam¬ 
ple,  the  effects  of  base  bleed  [see,  for  example,  [Sahu  et  al.  (1985)],  [Sahu  Sz  Heavey  (1995)],  and 
[Danberg  k  Nietubicz  (1992)]]  and  base  burning  [see,  for  example,  [Nietubicz  k  Sturek  (1988)]]. 
In  all  simulations,  the  results  show  a  correct  trend  for  the  drag  reduction  when  compared  with 
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experimental  results.  However,  in  some  cases,  even  global  flow  fleld  characteristics  such  as  the 
centerline  velocity  are  not  in  good  agreement  with  experiments.  This  indicates  that  there  are 
still  important  aspects  of  the  physics  which  are  not  captured  by  the  used  models.  Therefore 
empirically  determined  constants  have  to  be  introduced  which  need  to  be  adjusted  for  the  given 
geometry  and  thus  do  not  allow  prediction  of  flow  flelds  for  diflFerent  geometries. 

1.1.2  Investigations  of  the  Time-Dependent  Flow  Field 

All  investigations  mentioned  above  focused  on  calculating  the  time-averaged  turbulent  mean  flow 
in  the  near  wake  of  an  axisymmetric  body.  If,  however,  the  supersonic  waJce  flow  fleld  is  domi¬ 
nated  by  large  structures  (like  its  incompressible  counterpart),  it  is  necessary  to  investigate  the 
time  dependent  evolution  of  the  flow  field.  For  example,  the  time  dependent  behavior  of  a  two 
dimensional  supersonic  transitional  wake  has  been  investigated  by  [Chen  et  al.  (1989)]  [see  also 
[Chen  et  al.  (1990)]]  using  temporal  Linear  Stability  Theory  (LST)  and  Direct  Numerical  Sim¬ 
ulations  (DNS).  In  their  research  they  found  that  the  wake  is  subject  to  temporal  instabilities. 
Disturbances  in  the  shear  layer  get  amplified  and  form  large  structures  which  dominate  the  flow 
field  further  downstream.  Increasing  the  Mach  number  has  a  damping  efiect  on  these  instabili¬ 
ties.  It  is  as  yet  unclear  if  the  supersonic  axisymmetric  wake  flow  is  also  subject  to  instabilities 
(temporal  or  spatial,  convective  or  absolute)  similar  to  its  incompressible  counterpart. 

It  is  well  known  that  for  incompressible  wakes  the  dynamics  of  large  (coherent)  struc¬ 
tures  play  a  dominant  role  in  the  local  and  global  behavior  of  the  flow.  This  evidence  was 
found  from  both  experimental  investigations  and  numerical  simulations  and  was  confirmed  by 
theoretical  studies.  For  incompressible  flow  past  bluff  bodies,  it  has  been  well  established 
that  the  existence  of  absolute  and  global  instabilities  is  responsible  for  the  development  of 
the  large  structures  [see,  for  example,  [Huerre  &  Monkewitz  (1990)]].  Using  numerical  simu¬ 
lations,  absolute  and  global  instabilities  were  found  for  a  two-dimensional  bluff  body  with  a 
blunt  base  by  [Hannemann  &  Oertel  (1989)]  and  for  an  axisymmetric  body  with  a  blunt  base 
by  [Schwarz  (1996)].  The  absolutely  or  globally  imstable  modes  for  the  axisymmetric  case  are 
of  a  helical  nature.  For  compressible  wahes,  especially  at  supersonic  speeds,  however,  relatively 
little  is  known  about  the  dynamical  behavior  of  the  large  structures  in  turbulent  flows  or,  in 
paxticular,  if  absolute  or  global  instabilities  exist.  This  is  true  for  supersonic  flows  in  general 
and  for  axisymmetric  wakes  in  particular. 

Due  to  the  lack  of  guidance  from  experimental  investigations,  no  successful  computational 
or  theoretical  attempts  have  been  made  in  the  past  to  study  the  unsteady,  dynamical  behavior 
of  laminar,  transitional,  or  turbulent  supersonic  axisymmetric  wake  flows.  For  the  subsonic 
(incompressible)  case,  there  is  considerable  evidence  that  the  evolution  of  the  large  structures  is 
due  to  the  hydrodynamic  instability  of  the  (time-averaged)  mean  flow  and  that  the  development 
of  these  structures  can  be  captured  by  stability  theory.  Experimental  results  for  incompress¬ 
ible  turbulent  mixing  layers,  two-dimensional  turbulent  wakes,  and  axisymmetric  wakes  and 
comparisons  with  stability  theory  have  shown  that  certain  key  features,  such  as  dominant  fre- 
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quencies,  mode  shapes  (amplitude  distributions),  and  streamwise  spacing  (wave  lengths)  of  the 
structures  can  be  well  predicted  by  linear  stability  theory  [see,  for  example,  [Wygnanski  (1994)] 
and  [Marasli  et  al.  (1989)]].  These  investigations  support  the  notion  that  hydrodynamic  insta¬ 
bilities  give  rise  to  the  generation  and  development  of  the  large  structures.  The  dynamical 
behavior  of  these  structures  is  responsible  for  the  strong  unsteady  flow  behavior  in  the  wake. 
Thus,  in  reality,  there  is  no  steady  turbulent  wake  flow,  even  when  small-scale  (high-frequency) 
fluctuations  are  not  talcen  into  consideration.  Rather,  the  flow  is  highly  unsteady,  dominated 
by  large-amplitude  and,  relative  to  the  small  scales,  low-frequency  fluctuations.  Because  of  the 
nonlinear  interaction  between  the  various  fluctuation  components,  the  actual  mean  flow  may  be 
strongly  dependent  on  the  composition  of  the  fluctuating  parts  of  the  flow  fleld.  In  addition,  it 
has  been  shown  that,  with  artificial  forcing,  existing  components  (structures)  can  be  modified  or 
new  components  (structures)  created  and,  as  a  consequence,  that  the  mean  flow  can  be  strongly 
altered  [see,  for  example,  [Wygnanski  (1994)]  and  [Marasli  et  al.  (1989)]]. 

Due  to  the  difficulties  of  the  experiments,  until  recently  it  was  unclear  if  large  structures  are 
present  in  supersonic  wake  flows  or  if  they  play  a  similarly  dominating  role.  Recent  experimental 
findings  indicate  that  large  structures  indeed  exist  for  the  case  of  two-dimensional  supersonic 
wakes  [see,  for  example.  Smith  Dutton  (1968)].  Some  quantitative  evidence  of  the  existence 
of  dominant  large  structmres  in  supersonic  axisymmetric  wake  flows  has  been  provided  by  the 
experiments  of  Demetriades  (1968),  who  investigated  the  unsteady  nature  of  the  flow  field.  The 
amplitude  spectra  in  his  results  display  distinct  peaks  at  certain  (relatively  low)  frequencies, 
thus  indicating  the  presence  of  dominant  structures.  If  large  structures  are  indeed  dominat¬ 
ing  the  near  wake  flow  field,  it  is  not  smprising  that  the  mean  flow  fields  strongly  depend  on 
experimental  conditions  and,  even  for  the  same  facility,  can  vary  with  minor  changes  of  the  ex¬ 
perimental  parameters.  Also,  for  this  reason,  it  is  not  surprising  that  current  turbulence  models 
for  calculating  mean  flows  using  Reynolds-averaged  Navier-Stokes  formulations  are  performing 
poorly  for  these  flows.  Time  dependent  methods  like  Large-Eddy  Simulations  (LES)  would  be 
more  adequate  to  capture  the  strongly  time  dependent  behavior  of  the  flow  field.  Chances  to 
arrive  at  better  turbulence  models  (for  Reynolds-Averaged  Navier-Stokes)  are  rather  slim  un¬ 
less  the  physical  and  dynamical  behavior  of  the  large  structiures  are  better  understood  and  this 
knowledge  is  implemented  in  the  turbulence  modeling. 

1.2  Present  Research 

One  way  to  achieve  more  insight  into  the  dynamical  behavior  of  the  supersonic  axisymmetric 
wake  flow  is  to  investigate  the  development  of  disturbances  in  laminar,  transitional  and  turbulent 
flows.  By  determining  the  mechanisms  that  govern  the  behavior  of  large  disturbances  at  the 
onset  to  turbulence  one  might  gain  more  understanding  about  the  behavior  of  the  fully  turbulent 
flow.  This  can  be  further  supported  by  comparing  the  time  dependent  behavior  of  the  laminar 
flow  with  that  of  the  turbulent  flow.  As  a  first  step,  the  question  needs  to  be  answered,  if  an 
absolute  or  a  global  instability  also  exists  for  the  supersonic  wake  flow.  For  this  reason,  in  the 
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present  research  the  evolution  of  artificially  generated  disturbances  has  been  investigated  for 
laminar  and  for  turbulent  fiow  fields.  Direct  Numerical  Simulations  (DNS)  and  Large-Eddy 
Simulations  (LES)  have  been  employed  for  the  investigation  of  the  time  dependent  behavior  of 
the  supersonic  axisymmetric  wake  flow.  All  simulations  are  based  on  solving  the  complete  set  of 
the  compressible  Navier-Stokes  equations  (spatially  filtered  for  LES)  in  a  cylindrical  coordinate 
system. 


Figure  1.3:  Region  of  interest. 

The  region  of  interest  for  the  given  problem  is  the  near  wake  of  the  axisymmetric  body 
as  shown  in  figure  3.1.  The  last  part  of  the  boundary  layer  approaching  the  base  corner  was 
included  in  the  calculation.  This  allowed  the  introduction  of  controlled  disturbances  through 
a  circular  blowing  and  suction  slot  on  the  body  into  the  boundary  layer.  As  the  initial  step, 
a  laminar  axisymmetric  steady  flow  field  was  calculated.  This  flow  field  was  then  artificially 
disturbed  either  by  blowing  and  suction  through  a  circular  slot  on  the  body,  or  by  a  local  pulse 
within  the  recirculating  region.  Then  the  time  dependent  evolution  of  the  entire  flow  in  response 
to  the  artificially  generated  distmbances  was  monitored.  For  the  given  flow  situation,  air  can 
be  treated  as  a  calorically  perfect  gas,  resulting  in  a  constant  heat  capacity  (cp)  and  a  constant 
Prandtl  number  (Pr). 

One  of  the  main  interests  in  this  research  was  the  interaction  between  the  recirculation 
region  and  the  external  flow  which  are  connected  by  a  shear  layer.  As  mentioned  above,  it  is 
known  that  for  subsonic  incompressible  flows,  leurge  structures  form  in  the  shear  layer  and  play 
a  dominant  role  in  the  local  and  global  behavior  of  the  wake  flow  and,  consequently,  have  a 
strong  effect  on  the  drag.  However,  for  supersonic  compressible  wake  flows  the  role  of  large 
structmes  in  the  dynamical  flow  behavior  and  in  particular  their  influence  on  the  base  pressure 
is  not  well  understood.  The  ultimate  goal  of  the  present  research  is  to  contribute  towards  the 
understanding  of  the  dynamical  behavior  of  the  large  structures  and  their  effect  on  the  global 
flow  behavior. 


Chapter  2 


Governing  Equations,  Initial  and 
Boundary  Conditions 


As  pointed  out  in  the  literature  review,  previous  investigations  have  focused  on  finding  a  time- 
averaged  solution  for  turbulent  axisymmetric  wake  fiows  at  supersonic  speeds.  In  almost  all  these 
research  efforts  the  resulting  flow  field  was  assumed  to  be  axisymmetric.  Therefore,  typically 
a  set  of  the  compressible  axisymmetric  Reynolds-averaged  Navier-Stokes  equations  was  used 
as  the  governing  equations.  In  some  cases,  these  were  further  simplified  by  using  thin-layer 
approximations.  In  the  present  work  the  main  emphasis  was  on  the  time-dependent  evolution 
of  the  near  wake  flow  field,  which  does  not  necessarily  have  to  be  axisymmetric.  In  fact,  for 
the  case  of  the  subsonic  wake  flow  the  instantaneous  flow  field  and  sometimes  even  the  time- 
averaged  flow  field  does  not  exhibit  axisymmetry  [see,  for  example,  [Schwarz  (1996)]].  Therefore, 
the  governing  equations  chosen  for  the  present  investigation  are  the  complete  three-dimensional 
unsteady  compressible  Navier-Stokes  equatioirs.  Since  the  given  geometry  is  axisymmetric  the 
equations  are  solved  in  a  cylindrical  coordinate  system. 

Three  different  sets  of  the  equations  are  being  used  according  to  the  different  tasks.  For  the 
determination  of  laminar  axisymmetric  flow  fields,  the  equations  are  reduced  to  the  axisymmet¬ 
ric  compressible  Navier-Stokes  equations  (see  section  2.2).  The  fully  three-dimensional  flow  field 
for  lower  Reynolds  number  flows  (laminar  or  transitional)  is  calculated  by  solving  the  complete 
three-dimensional  unsteady  compressible  Navier-Stokes  equations  (see  section  2.3).  For  the  so¬ 
lution  of  turbulent  flows  at  higher  Reynolds  numbers  Large-Eddy  Simulation  (LES)  is  employed, 
where  the  three-dimensional  equations  are  filtered  in  space  in  order  to  resolve  only  structures 
that  are  larger  than  a  specified  threshold  (see  section  2.4). 


11 
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2.1  Nondimensionalization 

All  variables  in  the  given  problem  are  nondimensionalized.  The  length  scale  chosen  for  the 
problem  is  the  diameter  D  of  the  axisymmetric  body.  Thus, 


^  D  '  ^  D' 

The  velocity  components  are  nondimensionalized  by  the  free  stream  velocity  Uo 

V*x  V* 

=  =  ;  ^e  =  7j-- 

Uqq  ^OO  ^00 


(2.1a) 


(2.1b) 


Temperature,  density,  viscosity,  and  heat  conductivity  are  nondimensionalized  by  their  values 
in  the  free  stream: 


T  =  ;  p  =  —  ;  ^ 

■Loo  Poo  Poo  ^OO 


(2.1c) 


This  automatically  leads  to  the  following  nondimensionalization  for  time,  pressure  and  energy: 


-  P  _  ® 


(2.1d) 


The  global  Reynolds  number,  the  Prandtl  number  and  the  free  stream  Mach  number  are  defined 

by 

j^^^PccUooD  .  .  Moo  =  — ,  (2.1e) 


ifeg  =  '  ;  pr  =  !::^ 

/ioo  ^oo 

where  the  speed  of  sound  in  the  free  stream  is  given  by 

0,00  —  ■\/'yi?jToo, 


(2.1f) 


7=—  ;  Ri  =  Cp-Cv. 


(2.1g) 


The  coefficient  of  the  conduction  term  in  the  energy  equation  is  given  by  the  product  of  the 
Peclet  number  and  the  Eckert  number 


N  =  P^Ec  =  (7  -  VjPrl^DMl^, 


(2.1h) 


where 


Pe  =  RepPh-  ;  Ec  = 


2.2  Axisymmetric  Calculations 


2.2.1  Equations 

For  the  calculation  of  axisymmetric  flow  fields,  the  unsteady  axisymmetric  compressible  Navier- 
Stokes  equations  are  solved.  The  equations  are  used  in  conservative  formulation  and  are  ex¬ 
pressed  in  a  cylindrical  coordinate  system.  Then,  they  can  be  written  in  the  following  vector 
form 
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d  f*  d  ^ 


where 


A  = 


B  = 


/pN 

pVx 

pVr 

\pej 

PVx 


pVx+P  +  Txx 

P^X^T  "1“  '^TX 

er 


(2.2a) 


(2.2b) 


\Vx{pe  +  p)  -  A  VxTxx  +  VrTrx  ) 


(2.2c) 


pVr 

P^x^r  "b  '^rx 
pVr+P  +  Trr 

tj  ar 


y  Vr(pe  +  p)  N  dr  '^x'^rx  “I"  V^T^r  J 


(2.2d) 


and 


/ 


E  = 


pVr 

pVxVf  “h  ^TX 

pVr  +  ^rr  “  Tee 


(2.2e) 


\Vr{pe  +  p)  -  +  VxTrx  +  VrTrr  J 

The  stress  tensor  components  (ry)  for  axisymmetric  flow  are  defined  by  [see,  for  example, 
[Bird  et  al.  (I960)]] 


2  /i 


-  1  J.  J- 

3  Red  \  dx  dr  r 
—  ^  p  r^g  _  ^ 

3  Ebd  .dx  dr  r 


2  p 

~  3^ 


dvx  2^ 

9r  r 


and 


'^xr  —  Ttie  — 


dVx  dVr 
Rbd  \.dr  ^  dx 


(2.3a) 

(2.3b) 

(2.3c) 

(2.3d) 


The  energy  equation  is  used  in  the  form  of  a  conservation  equation  for  the  total  energy  (e), 
which  is  defined  by 

f  T. .  1 

(2.3e) 


CvT(x>  1  _  _ 


The  viscosity  of  air  is  assumed  to  be  a  function  of  temperature,  according  to  Sutherland’s  law 
[see,  for  example,  [White  (1991)]],  which  is  given  by 


p{T)  =  tI  .where  Soo  = 


lllK 


(2.4) 
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The  same  relationship  holds  for  the  temperature  dependence  of  the  thermal  conductivity,  because 
the  Prandtl  number  is  assumed  to  be  constant.  Since  air  is  treated  as  a  calorically  perfect  gas, 
the  equation  of  state  is  given  by  [see,  for  example,  [Thumm  (1991)]] 


PT 

JMI,- 


(2.5) 


2.2.2  Integration  Domain 

The  domain  for  numerical  integration  of  the  axisymmetric  equations  is  shown  in  figure  2.1.  It 
extends  from  a  point  xq  on  the  body,  past  the  base  to  a  point  xi  sufficiently  far  downstream  of 
the  base,  such  that  the  recirculation  region  is  fully  enclosed  (typically  3  to  20  diameters  behind 
the  base).  In  the  radial  direction  it  starts  at  the  axis  of  symmetry  (r  =  0)  and  extends  typically 
2  to  10  body  diameters  away  fi:om  the  body. 


Figure  2.1:  Integration  domain  for  the  axisymmetric  equations. 


2.2.3  Boundary  Conditions 

The  integration  domain  for  the  axisymmetric  problem  includes  five  physical  boundaries,  as  is 
shown  in  figure  2.1.  These  boundaries  and  the  boundary  conditions  that  are  prescribed  are 
listed  below. 

Inflow  boundary 

At  the  inflow  boundary  all  flow  variables  {p,  Vx,  Vr,  T,  and  p)  are  prescribed.  In  the 
subsonic  layer  the  pressure  (p)  is  extrapolated  from  the  field  assuming  a  constant  pressure 
in  streamwise  direction. 


Outflow  boundary 

At  the  outflow  boundary  the  second  derivative  in  axial  direction  for  all  primary  flow 
variables  (p,  pvx,  pvr,  and  pe)  is  set  to  zero.  This  boundary  condition  has  been  introduced 


CHAPTER  2.  GOVERNING  EQUATIONS,  INITIAL  AND  BOUNDARY  CONDITIONS  15 


by  Fasel(1976)  and  was  found  to  allow  disturbances  to  travel  through  the  boundary  without 
major  reflections. 

Walls 

At  the  walls,  all  velocities  are  set  to  zero. 


Ux  =  0  ;  Ur  =  0 


(2.6a) 


For  the  calculation  of  a  steady  flow  fleld  the  body  is  assumed  to  be  adiabatic.  Thus,  at 
the  horizontal  wall  it  follows 

f  =  0,  (2.6b) 

and  at  the  vertical  wall 

or 

&=»• 

For  the  calculation  of  an  unsteady  flow  field,  the  body  temperature  is  kept  constant.  Thus 


(2.6d) 


Axis  of  symmetry 

Since  the  flow  field  is  assumed  to  be  axisymmetric,  at  the  axis  of  symmetry  it  follows  that 


dvx  dp  de 

^=0  ;  ..  =  0  ;  ^  =  0  ;  ^=0. 


(2.6e) 


Free  stream  boundary 

At  the  free  stream  boundary,  three  difierent  boundary  conditions  are  applied,  depending 
on  the  free  stream  Mach  number  of  the  problem.  A  detailed  description  of  the  free  stream 
boundary  conditions  is  given  in  appendix  B. 


Initial  condition 

As  an  initial  condition  for  the  calculation  of  the  steady  axisymmetric  flow  field,  the  sim¬ 
ilarity  solution  for  a  compressible  flat  plate  boundary  layer  is  used.  The  radial  velocity 
distribution  is  calculated  by  integrating  the  steady  continuity  equation.  For  the  calculation 
of  an  unsteady  axisymmetric  flow  field,  a  previously  calculated  steady  flow  field  is  used  as 
initial  condition. 


2.3  Three-dimensional  Calculations 


2.3.1  Equations 

The  fully  three-dimensional  flow  field  is  calculated  by  solving  the  complete  unsteady  three- 
dimensional  compressible  Navier-Stokes  equations.  As  for  the  axisymmetric  flow  field,  the  equa¬ 
tions  are  expressed  in  conservative  formulation  in  a  cylindrical  coordinate  system.  Again  they 
can  be  written  in  the  vector  form 

^A  +  -^B  +  ^C+  --^D  +  -E  =  0, 
dt  dx  dr  r  do  r 


(2.7a) 
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where  now 


(2.7b) 


pvl+p  +  Txx 

pVxVr  +  Trx 

pVxVe  4-  TQx 

^Vx{pG  A  p)  ~  ^x'^xx  "b  ’^tT'tx  "b  '^B'^Ox  , 

'  pVr  '' 

pVx^r  “b  '^rx 
pV^+p  +  Trr 
pVrVe  +  TrB 

^Vr{pe+p)- VxTrx  +  VrTrr  +  V0Tre  j 


(2.7c) 


(2.7d) 


pvxVe  +  Tex 
pVrVe  +  TrB 
pVg+p  +  Toe 

,Ve{pe+p)  -  +  VxTBx  +  VrTre  +  VgTgg 


{2.7e) 


pVxVr  +  Trx 
pvl-pvl  A  Trr  -  Tgg 


(2.7f) 


2pVrVg  +  2Trg 

\vr{pe  +  p)  -  +  VxTrx  +  VrTrr  +  VgTrB  ) 

For  a  three-dimensional  flow  field,  the  stress  tensor  components  in  a  cylindrical  coordinate 
system  are  defined  as  follows. 

_  _  P  \^dvx  2,^  1  2  p  \  ^dvx  dvr  Idvg  Vr]  . 


__  p  \^dvr  _  2  1  _  2  _  5^  1  ^  V 

Ren  .dr  3  ^  .  Sifep  _dx  dr  r  86  r.  ’ 


(2.8b) 


Tee  = 


P  l^fldvg  Vr\  2  ]  2  P  \dVx  dVr  2  dvg  ^Vr]  . 


Txr  —  Trx  — 


p  \dvx  .  dv, 


Reo  \.dr  dx  \  ' 


(2.8d) 


CHAPTER  2.  GOVERNING  EQUATIONS,  INITIAL  AND  BOUNDARY  CONDITIONS  17 


— 


Rei)  [  dx 


dve  1  dvx 

r\  I  €-\/\ 


(2.8e) 


TQr  —  — 


d  \ve 


n  r  dvg  Vg  1  dVr 

Rej)  _dr  r  ^  r  86 


•  (2-8f) 


As  before,  air  is  treated  as  a  calorically  perfect  gas  and  the  total  energy  term  and  the 
temperature  dependence  of  the  viscosity  are  defined  in  the  same  way,  as  shown  in  equation  2.5 
(see  section  2.2). 


2.3.2  Integration  Domain 

The  integration  domain,  as  shown  in  figure  2.2,  includes  part  of  the  boundary  layer  on  the 
axisymmetric  body  and  extends  downstream  beyond  the  recirculation  region  (typically  3  to  20 
diameters  downstream  of  the  base). 


Figure  2.2:  Integration  domain  for  the  three-dimensional  equations. 


2.3.3  Boundary  Conditions 

As  can  be  seen  in  figure  2.2,  the  integration  domain  for  the  three-dimensional  problem  includes 
four  boundaries.  The  boundaries  and  the  prescribed  boundary  conditions  axe  listed  as  follows. 

Inflow  boundary 

For  the  three-dimensional  (unsteady)  calculations,  all  variables  at  the  inflow  boundary 
axe  specified.  The  only  exception  is  the  pressure  in  the  subsonic  layer,  where  the  first 
derivative  of  the  pressure  in  axial  direction  is  set  to  zero. 

Outflow  boundary 

As  for  the  axisymmetric  calculations,  at  the  outflow  boundary  the  second  derivatives  in 
axial  direction  of  all  primary  flow  variables  are  set  to  zero. 
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Walls 

At  the  walls  all  velocities  are  set  to  zero 

=  0  ;  Vr  =  0  ;  ve  =  0  (2.9a) 

The  temperature  on  the  walls  is  held  constant  based  on  the  assumption  that  the  body 
cannot  change  temperature  on  the  time  scale  of  the  disturbances.  Thus 

f  =  0  (2.9b) 

The  zero-velocity  boundary  condition  is  applied  at  the  wall  except  on  the  circular  section 
of  blowing  and  suction,  which  is  used  for  disturbance  generation.  This  is  described  in  more 
detail  in  section  3.7. 

Free  stream  boundary 

The  boundary  conditions  at  the  free  stream  are  similar  to  the  axisymmetric  case.  As  before, 
three  different  kinds  of  boundary  conditions  are  applied  at  the  free  stream  boundary, 
according  to  the  free  stream  Mach  number  (Mqo).  The  three  different  boundary  conditions 
are  described  in  more  detail  in  appendix  B. 

Initial  condition 

A  previously  calculated  steady  axisymmetric  flow  field  is  used  as  the  initial  condition  . 

2.4  Large-Eddy  Simulations 

2.4.1  Spatially  Filtered  Three-Dimensional  Equations  for  LES 

At  sufficiently  high  Reynolds  numbers  the  flow  field  eventually  becomes  turbulent.  In  that  case. 
Direct  Numerical  Simulations  become  very  memory  and  CPU  time  intensive.  In  order  to  be  able 
to  calculate  the  flow  field  for  a  turbulent  Reynolds  number  the  equations  therefore  cannot  be 
solved  directly,  because  of  the  lack  of  spatial  resolution.  In  most  applications,  however,  the  small 
tmbulent  scales  are  not  of  major  interest.  For  that  reason,  methods  have  been  developed  to 
model  these  small  scales.  In  almost  all  previous  research  efforts,  the  turbulent  compressible  wake 
flow  field  was  calculated  by  using  a  Reynolds-averaged  form  of  the  Navier-Stokes  equations  (see 
section  2.3).  These  calculations  always  lead  to  a  time-averaged  flow  field.  In  the  present  research, 
however,  the  main  focus  is  on  the  time  dependent  behavior  of  the  flow  field  and  the  evolution 
and  dynamics  of  large  structures.  Therefore,  in  the  present  research  the  method  of  Large-Eddy 
Simulations  (LES)  was  chosen  for  the  calculation  of  the  time  dependent  turbulent  flow  field. 
LES  has  been  successfully  applied  to  isotropic  turbulence  and  to  wall  bounded  flows,  such  as 
flat  plate  boundary  layers  and  duct  flow.  For  LES,  the  Navier-Stokes  equations  are  filtered 
in  space,  using  the  assumption  that  the  small  turbulent  scales  are  separable  from  the  larger 
structures  (eddies).  The  latter  can  then  be  calculated,  while  the  former  have  to  be  modeled.  In 
the  following  sections  the  equations  for  the  resolved  quantities  are  derived.  This  derivation  was 


CHAPTER  2.  GOVERNING  EQUATIONS,  INITIAL  AND  BOUNDARY  CONDITIONS  19 


done  for  the  most  part  by  [Israel  (1996)]  in  close  collaboration  with  [Speziale  (1996)].  Here,  the 
derived  equations  are  furthermore  adapted  to  the  axisymmetric  geometry. 

2.4.2  Spatial  Filtering 

For  the  filtering  of  the  Navier-Stokes  equations,  any  spatial  filter  of  the  form  G{x\  A)  with  filter 
width  A  can  be  used,  provided  that  it  has  the  following  properties: 


G(-i-oo),  G( — oo)  — ^  0  ; 

/+00 

G(Od^  =  1. 

■00 

The  filtered  value  of  property  /,  denoted  as  /,  is  defined  by  the  convolution 

/+0O 

G{x-C-,A)fiOd^. 

-OO 


(2.10) 

(2.11) 


(2.12) 


This  allows  the  decomposition  of  any  quantity  into  a  filtered  and  a  fiuctuating  part  by 

/  =  7+/'. 

For  spatial  derivatives  of  filtered  quantities  Ghosal  &  Moin  (1993),  for  example,  pointed  out 
that 


dfjx) 

dx 


d  r+°° 

=  I:  /  Gix-  onm 

dx  J —00 

/+O0 

G'ix-OfiOd^ 

-00 

.  r+oo 

=  -  0(x  -  „  +  /  G(i  - 

J—oo 


^  df{x) 
dx  ’ 

since  /  must  be  finite  on  the  boundaries  and  the  filter  vanishes  at  infinity.  This  property  is  only 
valid  for  filters  with  a  uniform  filter  width.  For  non-uniform  filters,  it  can  be  shown  that  this  is 
only  true  to  a  second  order  approximation  [see,  for  example,  [Ghosal  &  Moin  (1993)]]. 

In  the  conservative  formulation  of  the  compressible  Navier-Stokes  equations,  the  velocity 
always  occurs  in  products  with  the  density.  When  the  equations  are  filtered  this  leads  to  the 
velocity-density  correlation  terms,  pui.  In  order  to  obtain  separate  variables  (that  is,  not  under 
a  correlation),  the  so  called  Favre  filter  (denoted  by  a  tilde)  is  introduced,  which  is  defined  by 


(2.13) 


The  Favre  filtered  variable  is  also  sometimes  referred  to  as  a  mass  averaged  quantity.  This  allows 
the  separation  of  the  density- velocity  correlation  terms  by  means  of  the  substitution  pf  =  pf. 
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(Note:  In  the  incompressible  limit,  p  becomes  a  constant  and  /  /.)  As  before,  the  physical 

quantity  can  also  be  separated  into  a  Favre  filtered  part  and  a  fluctuating  part  (denoted  by  a 
double  prime), 

/  =  /  +  /".  (2.14) 


2.4.3  The  Continuity  Equation 

The  continuity  equation  for  a  compressible  fluid  can  be  written  in  index  notation  as 

dp  dpuj  __ 
dt  dxi 

Applying  the  spatial  filter  to  this  equation  directly  leads  to  the  filtered  continuity  equation, 

dp  ,  dpu; 


(2.15) 


+ 


=  0. 


dt  dxi 

Introducing  the  Favre  filtering  the  correlation  in  the  second  term  can  be  eliminated  and  the 
continuity  equation  for  the  mass  averaged  velocity  can  be  obtained  as 

dp  dpui 


dt  ^  dxi 


=  0. 


(2.16) 


2.4.4  The  Momentum  Equations 

Similar  to  the  continuity  equation,  the  momentum  equations  for  a  compressible  fluid  can  be 
written  in  index  notation  as 


dpui  dpuiUj  _ 
dt  dxj 


dp 

dxi 


dxi 


ij 


dxj  ’ 


(2.17) 


where  the  stress  tensor  is  defined  by 


Tij  — 


p  I 


2duk 

Sdxr\ 


Re  ydxj  dxi 

Spatial  filtering  as  for  the  continuity  equation  leads  directly  to  the  filtered  momentum  equation 


dpui  ^  dpuiUj  _ 


dt 


dx  -i 


dp 

dxi 


dxj' 


(2.18) 


where  the  terms  UiUj  and  Tij  are  not  closed  and  have  to  be  modeled. 

To  preserve  the  form  of  the  unfiltered  equations,  the  term  uiuj  needs  to  be  replaced  with 
UiUj.  First,  the  term  can  be  expanded  as  follows  by  substituting  2.14, 


UiUj  =  {ui  +  u'l){uj  +  u'j) 

=  UiUj  +  Uiu'j  +  u'luj  +  u'lu'j, 


UiUj  =  UiUj  +  UiU'j  +  u'luj  +  u'iU'j. 


or  for  the  filtered  term 
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All  the  terms  containing  unknown  correlations  can  be  moved  to  the  right  hand  side.  To  combine 
them,  a  subgrid-scale  stress  tensor  is  defined  as  follows 


(Tij 


=  Lij  +  Cij  +  Rij 


(2.19) 


where  Lij,  Cij,  and  Rij  are  typically  called  the  Leonard  stress,  the  cross-stress,  and  the  Reynolds 
stress,  respectively.  They  are  defined  by 


Lij  —  p^Uitlj  UiUj'j, 
^ij  =  —p{uiu"j  Au'luj), 
Rij  =  -pu'fu'!. 


Substituting  2.19  into  the  filtered  momentum  equation  2.18  leads  to  the  mass-averaged  (Favre 
filtered)  momentum  equation 


&pui  dpuiUj  _  dp  d  . , 

St  dxj  Sxi  dxj 


(2.20) 


The  Leonard  stress  can  be  calculated  explicitly.  The  cross-stress  and  Reynolds  stress  are  not 
closed,  and  must  be  modeled  (see  section  2.4.7). 

To  compute  the  filtered  viscous  stress  tensor,  rij  needs  to  be  expressed  in  terms  of  the 
filtered  quantities.  In  many  applications,  the  molecular  viscosity  is  assumed  not  to  vary  due 
to  the  turbulent  fluctuations.  For  example,  [Erlebacher  et  al.  (1992)]  and  [Zang  et  al.  (1992)] 
both  assume  small  temperature  fluctuations,  and  therefore  take  /z  to  be  constant  over  the  filter 
width.  This  leads  to  the  following  expression  for  the  stress  tensor 


—  =  At(^)  ( ^  ^  2  \ 

Rbd  Sxi  3  dxk  y  ’ 


(2.21) 


where  the  molecular  viscosity  fj,{f)  is  a  function  of  the  Favre  filtered  local  temperature  only. 


2.4.5  The  Energy  Equation 


In  the  formulation  of  the  Navier-Stokes  equations  used  in  the  present  research  the  energy  equa¬ 
tion  is  expressed  in  terms  of  total  energy.  The  equation  in  index  notation  is  given  by 


dpe 

dt 


dxi 


,  1?  dT 

-u,{fe+p)  +  -— 


UjTij  . 


(2.22) 


Unfortunately,  if  this  equation  is  filtered  directly  several  terms  cannot  be  calculated  directly  and 
need  to  be  modeled.  This  includes  the  filtered  total  energy,  defined  as 
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where 


"  t;2  • 

This  quantity  is  not  closed  since  u^i  is  not  known.  However,  one  can  .construct  a  quantity  Er, 
which  here  is  called  resolved  scale  energy,  given  by 


Er  =  p  {CvT  +  . 


(2.23) 


This  quantity  can  be  directly  calculated  from  the  Favre  filtered  velocities  and  temperature.  In 
addition,  it  can  be  related  to  Et  using  the  relation 

Et  =  Er  —  2^ii‘ 

The  governing  equation  for  Er  can  be  derived  in  the  same  manner  as  the  equation  for  Et-  First, 
the  equation  for  the  resolved  kinetic  energy  ttjUj  can  be  derived  by  multiplying  the  resolved  scale 
momentum  equation  2.20  with  Uj,  which  leads  to 


_  dpui  .  dpuiUj  _  ~  d  j _  , 


dpuiUi  dfpUiUiUj 
dt  ^  dxj 


Ui  p-—  +  fnij—  -u 


dp  .  d  , 


'dxi  ^dx 


Ui  n  ['^ij  <^ij]  ■ 


Using  the  continuity  equation,  2.16,  multiplied  by  it  follows 

dpujUj  ^  dpujUjUj  _  .  \dpui  dpuiuA  _  dp  .  d 


dt  dxi 


The  right  hand  side  can  now  be  rewritten  using  the  momentum  equation  2.20  to  yield  the 
equation  for  the  resolved  scale  turbulent  kinetic  energy 


1  dpUiUi  1  dpUiUiUj 

2  dt  2  dxj 


_  -  (  ^  dr—  A 

dxi 


(2.24) 


Similarly,  the  temperature  equation  for  a  compressible  fluid  in  index  notation  is  given  by 


dpCyT  dpC-oTu. 
dt  dxi 


i  dUi  d  (  d  dT\  dUj  m  nr\ 


[see,  for  example  [Bird  et  al.  (I960)]]. 

Again,  the  equation  can  be  filtered  in  a  similar  manner  as  the  momentum  equation.  This 
leads  to 


dpCyT  dpCyTui 


dui  d  (  d  dT 


=  -P^  + 


dxi  dxi  V  N  dx 


-  T, 


^^dx^- 


(2.26) 


Similar  to  the  momentum  equation,  the  term  Tu,  cannot  be  calculated  directly  and  has  to  be 
modeled.  Introducing  a  subgrid-scale  heat  flux,  corresponding  to  the  subgrid-scale  stress  tensor 
2.19,  one  can  define 

Qi  =  -pCvifui  -  Tui) 

=  Lf  +  c9  +  R9^ 
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where 

if  =  pCvifui  -  fui), 
C?  =pC,{f^'!-THn), 
and 

r9  =  pC„f^l 


are  defined  similarly  to  the  Leonard  stress,  the  cross  stress  and  the  Reynolds  stress  tensors, 
respectively.  With  these  definitions  the  filtered  energy  equation  in  terms  of  the  Favre  filtered 
temperature  becomes 


dfpC.oT  dpCyTui  _  dui  d  df 
dt  dxi  ^  dxi  dxi  N  dxi  * 

where  the  right  hand  side  is  still  not  closed  without  additional  assumptions  to  model  the  filtered 
correlations  (see  section  2.4.7). 

Now  the  equation  for  resolved  scale  energy,  Er,  can  be  determined  by  summing  2.24  and 
2.27,  and  rearranging  indices.  This  leads  to 


duj 


(2.27) 


SEr  duiER 


dt 


dxi 


=  -p 


dui 


dxi  dxi 


d  dT 
Ndxi 


Qi 


-  Ti 


dui  . 


_  Ar- 
y  dxi  dxj 


This  equation  cannot  be  written  in  a  form  similar  to  2.22  due  to  the  presence  of  disparate  types  of 
filtering  in  the  pressure  and  molecular  stress  terms.  However,  following  [Erlebacher  et  al.  (1992)] 
the  pressme  terms  can  be  rewritten  as 


dui  _  dp 
^  dxi  dxi 


dpui  f_dui  dui 
dxi  ^dxi^ 

dpu 


dxi 

dpu. 

dxi 


I  l-dui  _dui  du'l  ,dui  ^du'! 

r  axi  ^dxi  dxi^^  dxi'^'^  dxi 


—  € 


V 


In  addition,  according  to  [Erlebacher  et  al.  (1992)]  the  terms  combined  in  are  negligible. 
The  viscosity  term  can  be  handled  similarly, 


dui  _  dfij  _  duiTij  I  dui  _ dui\ 

dxj  dxj  dxj  y  dxj  dxj  j 

^  dUjTij  (  dUj  dUi  du'l  dUj 

dxj  y  dxj  dxj  dxj  dxj  dxj  J 

_  dUiTij  ^ 

—  —3  r  ^Ti 
dXj 

where  the  terms  combined  in  e,-  are  neglected  for  similar  reasons  as  for  Cp  [see,  for  example, 
[Erlebacher  et  al.  (1992)]]. 
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The  subgrid-scale  stress  term  can  be  written  as  (noting  that  aij  is  symmetric): 


_  da. 


Ui 


d  ,  duj 

ax, 


11  — 


9  .  if  duj  dui 

9.-1  X 
~  9xi  ~ 


where  Sij  is  the  strain  tensor  formed  with  the  Favre  filtered  velocities. 

Introducing  all  the  above  assumptions  and  definitions,  the  final  form  of  the  resolved  energy 
equation  then  becomes 


dpR 

dt 


dxi 


—Ui{Eji  +p)  + 


\Ndxij 


Qi  Uj^Tij  aij) 


aij  Sij . 


(2.28) 


2.4.6  Equation  of  State 


In  addition  to  the  Navier-Stokes  equations,  the  equation  of  state  also  needs  to  be  filtered. 
Recalling  the  equation  of  state  for  an  ideal  gas 


„  PT 
^  iMl' 

it  follows  that  the  filtered  equation  is 


(2.29) 


(2.30) 


2.4.7  Subgrid-Scale  Model 


As  an  initial  step,  the  subgrid  model  chosen  for  the  present  research  is  a  constant  coefficient 
Smagorinsky  type  model  as  suggested  by  [Speziale  (1996)]  [see  also  [Smagorinsky  (1963)]].  The 
Smagorinsky  model  is  an  eddy  viscosity  type  model,  as  are  most  turbulent  models  currently 
used  in  ffuid  dynamics  simulations.  In  this  initial  study  the  isotropic  part  of  the  turbulent  stress 
tensor  is  assumed  to  be  negligible.  Thus,  the  turbulent  stress  tensor  takes  the  form 

aij  =  2pi  {2SkiSki)^ {Sij  —  ’^S/ikSij),  (2.31) 

where  the  turbulent  length  scale  £  is  given  by 

where  Aj  is  the  filter  width  in  the  i-th  coordinate.  Here,  Cs  is  the  so  called  Smagorinsky  constant. 
At  this  point  the  most  appropriate  value  for  this  constant  for  the  given  fiow  problem  is  unclear. 
It  has  been  found  for  other  applications  that  this  constant  can  vary  between  approximately  0.05 
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and  0.2,  where  the  lower  values  have  been  used  for  wall  bounded  flows  and  the  higher  values 
have  been  used  for  isotropic  turbulence  [[Speziale  (1996)]]. 

In  addition,  for  wall  bounded  flows  a  wall  damping  function  was  suggested  by  [Speziale  (1996)], 
which  ramps  down  the  efiect  of  the  model  in  the  wall-near  region  (For  the  present  research  the 
turbulence  model  was  only  used  downstream  of  the  blunt  base.  Thus,  in  this  case  the  wall  of 
interest  is  the  base  itself.) 


1 

3 


where  x  is  the  wall  normal  direction  and 


A+  =  25  ;  =  ;  Ur=J^ 

p{T)  V  P 


x=0 


For  the  first  investigations,  the  constant  was  chosen  to  be  ^ 


Cs  =  0.065 

For  the  subgrid-scale  heat  stress,  a  gradient  transport  hypothesis  is  used  [see,  for  example 
[Erlebacher  et  al.  (1992)]],  which  leads  to  the  following  model 

=  (2.32) 

where  the  turbulent  Prandtl  number,  Prr  is  assumed  to  be  unity  for  the  present  research.  Using 
a  turbulent  conductivity  ‘dt,  defined  as 

=  ^^yl2SkiSki,  (2.33) 

the  subgrid  heat  stress  can  be  expressed  by 

Qi  =  (2.34) 


2.4.8  Resolved  Scale  Equations 


In  the  cylindrical  coordinate  system  that  has  been  used  here,  the  final  equation  for  the  resolved 
(filtered)  quantities  is  given  by 


^A+^B  +  ^C+-^D  +  -E  +  F=^0. 
at  ox  or  r  o9  r 


(2.35a) 


^Note:  =  Cr,  or  Cr  w  0.8409\/cr  where  Cr  is  the  Smagorinsky  coefficient  as  defined  in 

[Erlebacher  et  al.  (1992)]. 
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Here  the  vectors  are  defined  by 


pvl+P  +  Txx-  Oxx 


pVx^r  'A  Tfx  (Jfx 

p}xVe  +  Tgx  -  aex 

Vx{Er+p)  -  (I'  +  +  Vx{Txx-(^xx)  +  Wr(Tri-CTrx)  +  h{Tex-Oex) 


pVxVr  "I"  Trx  (^rx 
PV^  +p  +  Trr  -CTrr 
pUrVe  +  Tr9  -  (TrO 

MEr+p)  -  (^+  +  Vxirrx-Orx)  +  VrijTT-f^rr)  +  Ve{Tre-Ore) 


pVg 

pvxve  +  Tgx  —  aox 
pVrVg  +  Trg  -  ar9 
pVg  +P  +  Tgg  -  Ogg 

Vg{ER  +p)  -  +  Vxij9x-<yex)  +  Vr{rr9-(^r9)  +  M'^e9-(^9 


P^x'^T  '^rx  ^rx 
P^r~P^9+  ^rr  "  <^rr  "  T99  -  a9g 
2pVrVg  +  2Tr9  -  (7r9 

Vr{ER  +  p)  ~  +  ^x{Trx~(^rx)  d"  fir('7'rr “"O’rr)  d"  ^9{Tr9~^r9) 


pe^SW  [||5||"-|(Sxx  +  5rr  +  5,fl)' 


where  the  norm  of  the  strain  tensor  is  defined  as 


||S||  =  V2JP,,  +  S},  +  4  +  2S?.  +  2S?,  +  24, 


(2.35b) 


,  (2.35c) 


,  (2.35d) 


,  (2.35e) 


(2.35f) 


(2.35g) 


(2.36) 
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and  the  strain  tensor  components  are  given  by 


Q  dvx 

dx  ’ 

(2.37a) 

a  9Vr 

dr  ’ 

(2.37b) 

(2.37c) 

q  Q  1  ,  dVr' 

2  [dr  dx. 

,  (2.37d) 

Q  _  q  _  ^  Idvx 

2[dx^rde 

,  (2.37e) 

and 

q  _  c  1  L  ^  1  ^"’•1  ^\dve  ve  1  dvr]  ,  .. 

bgr  —  br6  —  — T  r—  —  H - —  —  — -  — - 1 - —  .  (2.37f) 

2  1  dr  [r  i  r  da  \  2  i  dr  r  r  dO  . 

Then  the  components  of  the  filtered  stress  tensor  can  be  written  as 

,  (2.38a) 

(2.38b) 

(2.38c) 

=  _ ■:=  _ n  q 

!  XT  —  *  TX  —  ^  rt 

I^D 

(2.38d) 

Tgx  —  “^xO  —  2  SxO, 

Ren 

(2.38e) 

and 

Tflr  —  TrO  —  ^  ^rOt 

(2.38f) 

and  the  components  of  the  tmbulent  stress  tensor  are  given  by 

CTxx  =  2pf\\S\\  (Sxx  + 

1  ,  (2.39a) 

arr  =  2pe^S\\  (^rr  +  ^V-^) 

(2.39b) 

age  =  2pf\\S\\  + 

1  ,  (2.39c) 
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||5||<Srx) 

(2.39d) 

(^6x  -  (^xB  =  2pI^\\S\\S:ce, 

(2.39e) 

and 

(^Or  =  (^re  =  2p^^||<S||5rfl, 

(2.39f) 

where  the  turbulent  length  scale  i  is  given  by 

£  =  CsA, 

(2.39g) 

and  the  filter  width  A  for  a  cylindrical  coordinate  system  by 

A  =  (Aa;ArrAe)5 . 

(2.39h) 

2.4.9  Integration  Domain  and  Boundary  Conditions 

The  integration  domain  and  the  boundary  conditions  for  the  Large-Eddy  Simulations  are  the 
same  as  for  the  three-dimensional  Direct  Numerical  Simulations  (see  sections  2.3.2  and  2.3.3). 


I 


Chapter  3 

Numerical  Method 


For  both  the  Direct  Numerical  Simulations  (DNS)  and  the  Large-Eddy  Simulations  (LES),  the 
governing  equations  and  the  boundary  conditions  axe  discretized  using  fourth-order  accurate 
difference  approximations  in  the  axial  (x)  and  radial  (r)  directions.  To  achieve  higher  resolution 
inside  the  shear  layer,  a  geometric  stretching  is  applied  in  both  directions.  Because  of  the  truly 
periodic  nature  of  the  flow  fleld  in  the  azimuthal  direction  a  pseudo-spectral  approximation  is 
employed  in  6,  using  a  truncated  Fourier  series.  According  to  the  pseudo-spectral  method,  the 
non-linear  terms  are  calculated  in  physical  space  and  then  transformed  into  spectral  space.  This 
method  is  based  on  the  numerical  method  used  by  [Thumm  (1991)].  Here,  the  discretization 
of  the  governing  equations  is  only  shown  for  the  three-dimensional  equations  of  the  Direct 
Numerical  Simulation.  The  axisymmetric  equations  are  just  a  subset  of  those  and  for  the  Large- 
Eddy  Simulations  the  equations  of  the  resolved  scales  are  very  similar. 

3.1  Discretization  of  the  Governing  Equations 

The  numerical  method  for  solving  the  governing  equations  consists  of  a  fourth-order  accurate 
explicit  Runge-Kutta  scheme  for  the  time  integration,  and  fourth-order  accurate  central  flnite 
differences  for  the  approximation  of  the  spatial  derivatives  in  the  radial  and  axial  directions. 

3.1.1  Fourier  Decomposition 

Due  to  the  periodic  nature  of  the  flow  fleld  a  truncated  Fourier  series  transform  is  applied  in  the 
circumferential  direction.  Thus,  the  flow  quantities  in  physical  space  can  be  expressed  by  their 
representation  in  Fourier  space  as  follows: 

f{x,r,e)=  ^k{x,r)  (3.1) 

k--K 
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where  i  =  y/^. 

Since  /  has  to  be  real,  it  follows  that 

F_k  =  conj.  .  (3.2) 

Substituting  the  Fomrier  representation  into  the  governing  equation  and  using  the  fact  that 


one  gets  a  set  of  +  1  coupled  equations  of  the  following  form 

Pi  ri  Pi  Pi  i  z! 

+  ^Ck  +  -Dk  +  -Ek  =  0  (3.4) 

ot  az  or  r  r 

The  coupling  terms  arising  from  cross  products  in  the  convective  and  the  dissipative  terms  are 
calculated  using  a  pseudo-spectral  method,  where  the  products  are  calculated  in  physical  space 
and  then  transformed  back  into  Fourier  space  (see,  for  example,  [Tourbier  (1991)]). 

The  Fomrier  representation  reduces  the  original  three-dimensional  problem  to  a  set  of  if  -I- 1 
two-dimensional  problems.  The  integration  domain  for  this  problem  set  is  the  same  as  for  the 
axisymmetric  problem,  shown  in  figure  2.1.  This  results  in  an  artificial  boundary  at  the  axis  of 
symmetry  that  was  not  present  in  the  original  three-dimensional  problem.  For  this  boundary, 
an  additional  set  of  boundary  conditions  has  to  be  specified.  This  will  be  further  explained  in 
section  5.2.2. 

3.1.2  Spatial  Derivatives 

In  order  to  enhance  numerical  stability,  the  fourth-order  accurate  central  finite  differences  that 
are  employed  for  the  approximation  of  the  spatial  derivatives  are  split  into  one-sided  backward 
and  forward  finite  differences.  The  splitting  directions  are  alternated  at  intermediate  time  steps. 
The  method  is  similar  to  the  MacCormack  method  [see,  for  example,  [MacCormack  (1969)]]. 
This  way  artificial  viscosity  is  added  to  the  discretized  equations,  damping  grid  mesh  oscillations. 
The  resulting  difierence  formulation  is  shown  the  following  example. 

The  fourth  order  central  difierences  are  split  into  a  second  order  backward  and  a  second 
order  forward  difference  in  the  following  manner: 

=  ;  //,  =  ^(-7/i  +  8/i+i-/i+2).  (3.5) 

Adding  both  derivatives  leads  to 

2/'  =  ^  (/i-2  -  8/i-i  +  8/i+i  -  fi+2) ,  (3.6) 

Dfl 

which  is  the  same  formula  as  the  commonly  used  fourth  order  central  difference,  except  that 
here  the  one  sided  differences  are  evaluated  at  different  intermediate  time  steps. 
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Second  derivatives  are  separated  into  two  first  order  derivatives,  with  the  splitting  applied 
in  the  opposite  direction  for  the  first  derivative  inside,  as  follows: 

fbi  =  ^  ~  +  ffi-2) 

=  ^  i-7fi-2  +  64/i_i  -  114/i  +  64/i+i  -  7fi+2) ,  (3.7) 

and 

ffi  =  ^  +  S/di+i  -  fbi+2) 

=  ^  (-7/i-2  +  64/i_i  -  114/i  +  64/hi  -  7/h2)  •  (3-8) 

After  adding,  this  leads  to  the  following  approximation  for  the  second  derivative 

2/"  =  i-7fi-2  +  m-i  -  114/i  +  64/hi  -  7/h2)  .  (3.9) 

Applying  Taylor  series  to  the  difference  formula  shows  that  the  approximation  is  fourth-order 
accurate: 

2/f  =  +  32/"(a;i)  32f"{xi)  -  Uf"{xi))  +  0{h^) 

=  2f"{xi)  +  Oih%  (3.10) 

Since  the  formulation  for  the  second  derivative  is  the  same  for  all  intermediate  time  steps,  no 
additional  artificial  viscosity  is  introduced. 

3.2  Grid  Stretching 

As  mentioned  above,  a  geometrical  stretching  is  applied  in  the  axial  and  radial  directions.  For 
the  stretching  functions  in  both  directions,  higher  order  polynomials  have  been  chosen  [see, 
for  example,  [Tourbier  (1991)]].  Through  the  grid  stretching,  extra  factors  are  added  to  the 
derivatives  of  the  original  equations.  Changing  from  the  {x,r)  coordinate  system  to  the  (C,r/) 
coordinate  system,  these  factors  are  given  by  ^  and  ^  for  the  axial  and  the  radial  direction, 
respectively.  The  stretching  functions  that  were  chosen  here  and  the  factors  are  given  in  the 
following  two  sections.  A  sample  of  the  stretched  grid  is  shown  in  figme  3.1. 

3.2.1  Axial  Direction 

For  the  grid  stretching  in  the  axial  direction,  a  third  order  polynomial  was  chosen  as  the  stretch¬ 
ing  function.  The  polynomial  is  of  the  following  form: 

x(0  =  +  Cx.  (3.11) 

To  guarantee  the  highest  resolution  to  be  at  the  base  (a;  =  0),  the  equation  has  to  satisfy  the 

condition  ^  =  0.  For  convenience,  the  following  conditions  were  added: 

“t  1=0 
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Figure  3.1:  Example  for  grid  stretching. 

ata:  =  0:  ;  ^  =  0 

at  a:  =  1:  ^  =  1. 

Here,  Sx  is  a  stretching  parameter  that  can  range  between  zero  and  one,  where  s®  =  0  means 
infinite  stretching  and  Sx  =  1  means  no  stretching.  Using  these  conditions,  the  three  constants 
for  the  stretching  function  are  given  by 

Ax  =  l  —  Sx  ;  Bx  =  Sx  ;  Cx  =  0.  (3.12) 

Thus,  the  stretching  in  the  axial  direction  is  given  by 

x(0  =  (1  -  Sx)^^  +  SxC 

The  additional  factor  that  is  added  to  the  axial  derivatives  is  given  by 

3(1-Sx)e  +  Sx 

3.2.2  Radial  Direction 

In  the  radial  direction  a  fifth  order  polynomial  was  chosen  for  the  stretching  function,  because 
the  point  of  highest  resolution  is  at  r=  0.5  and  the  function  needs  to  have  two  inflection  points. 
The  general  form  of  the  polynomial  is  given  by 

r(T})=ArTJ^ABrV^  +  Crr]  +  Dr.  (3.15) 

In  order  to  assure  the  highest  resolution  at  r  =  0.5,  the  polynomial  needs  to  satisfy  the  condition 

^  =  0.  For  convenience,  the  following  conditions  were  added: 

x=0.5 


(3.13) 

(3.14) 
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at  r  =  0: 
at  r  =  0.5: 


T]  =  0 

dr 

dr} 


1=0.5 


=  S. 


r  ) 


7?=  1. 


Again,  Sr  is  a  stretching  parameter  that  can  range  between  zero  and  one,  where  Sr  =  0  means 
infinite  stretching  and  Sr  =  0.5  means  no  stretching.  Using  these  conditions,  the  three  constants 
for  the  stretching  function  are  given  by 

and  the  stretching  function  in  radial  direction  is 

Hn)  =  -g  {s,  -  j)  »=  +  i  (*r  - 1)  +  [5  -  5  (“r  -  5)]  I-  (317) 


The  stretching  factor  for  the  radial  derivatives  is  given  by 


dr] 


iv)  = 


1 


dr  -f  (^Sr  -k)v^  +  T  («r  -  5)  +  [l  -  I  -  ^)] 

3.3  Discretized  Equations 


(3.18) 


The  resulting  discretized  equations  are  shown  here  for  the  example  of  the  continuity  equation 
for  the  time  t  =  to  + 1  At  at  the  point  ^  and  rj  =  r]o  +  nAr}.  The  continuity  equation 

for  the  A:— th  Fourier  mode  in  the  stretched  coordinate  system  is  given  by 


dh  ^  dj  d{fm)f.  _  d-q  d{pv)^.  _  ik{pw),^  _  {pv)t, 
dt  dx  dr  dq  r{q)  r{q) 

For  this  equation  the  discretized  Runge-Kutta  intermediate  steps  are: 
First  Predictor; 

Forward  Euler  half  step  with  backward  spatial  finite  difierences 

First  Corrector: 

Backward  Euler  half  step  with  forward  spatial  finite  difierences 

Pk  \m,n,l+\  ~  Mm,n,l  "I"  ~2  \m,n,l+\  ' 

Second  Predictor: 

Midpoint  rule  full  step  with  backward  spatial  finite  difierences 

RHS”L,„,,+J  ■ 

Second  Corrector: 

Simpson’s  rule  full  step  with  forward  spatial  finite  difierences 

Pk\m,n,l+l  ~  Pk\m,n,l  ^ 

+2  RHS"U,„,,+i  +  RHS"-U, 


(3.19) 


(3.20) 


(3.21) 


(3.22) 


(3.23) 


CHAPTER  3.  NUMERICAL  METHOD 


34 


where  the  ride  hand  sides  for  the  intermediate  time  steps  are  given  by 


RHS 


1 

6A^  dx 

■  ^ 

m 

7  {fnj)k 

—  8  {pu)u 

m,n,l 

+ 

Im— 2,Ti,i 

1  dr] 

- 

- 

6A77  dr 

n 

7  {pv)k 

-  8 

m,n,l  ” 

+  (P^)fc 

ik 

r{ri) 


m,n,l 


r(T)) 


(3.24) 
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RHS* 


lm,rt,i+5 
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6A^  dx 
1  dr] 


QArj  dr 
ik 


-7  Mfc 
-7  {fyv)l 


m,n,J+5 


1  +  8  {pu)k 


r{T]) 


ipw)k 


*1  mI 

m,n,J+5 


1+8  ipv)k 

m,n,i+5 


-  (pw)I 

771,71+1,^+5 


m+2,7i,i+^J 


m,n+2,i+| 


r(r?) 


(3.25) 


and 


RHS-  = 

1  ^ 
6A^  da: 

1  dr] 
6At]  dr 

ik 


7  (H* 


8  (pw)* 


77l“l,7l,/+5 


1  +  iP^)k 


m,n,i+5  m,n-l,Z+5 


m-2,n,/+5j 

m,n-2,/+i 


r(,?) 


(PW')*: 


(pw) 


m,n,i+“ 


^  l7n,n,i+“ 

r(?7) 


TJXJC***  I  _ 

lm,n,i+l  — 
1  ^ 
6A^  da; 

1  dr] 
6Ar]  dr 

ik 


-7  (p«)fc 

r-^v 

-7  (pv)fc 


m,n,l+l 


m,n,l+l 


+  8  (p«)jfe 

+  8  {pv)k 


m+l,n,l+l 


7n,n+l,J+l 


(pn) 

-  {pv)T 


*** 
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(3.26) 


m+2,n,/+lj 

m,n+2,i+l 


r{r]) 


{pw) 


*** 

k 


(pv) 


m,n,i+l 


^  1777,71, /+! 

r(77) 


(3.27) 


3.4  Integration  Domain 

For  the  discretized  governing  equations,  the  integration  domain  (shown  in  figure  3.2)  is  two- 
dimensional.  The  whole  domain  is  divided  into  two  subdomains.  The  first  part  includes  part  of 
the  boundary  layer  on  the  axisymmetric  body  and  extends  from  the  infiow  boundary  to  the  base 
corner  and  fi:om  the  wall  of  the  axisymmetric  body  to  the  free  stream  boundary.  Adjacent  to  the 
first  domain,  the  second  part  extends  from  the  base  to  a  location  downstream  of  the  recirculation 
region.  The  second  domain  extends  in  the  radial  direction  firom  the  axis  of  symmetry  to  the  free 
stream  botmdary.  For  the  three-dimensional  problem  this  adds  an  additional  boundary  at  the 
axis  and  boundary  conditions  have  to  be  specified  (see  the  following  section). 


3.5  Discretized  Boundary  Conditions 

As  in  the  previous  chapter,  the  boundary  conditions  for  the  discretized  problem  are  divided  into 
the  conditions  for  the  axisymmetric  case  and  the  conditions  for  the  three-dimensional  problem. 
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Figure  3.2:  Integration  domain  for  the  discretized  equations 
3.5.1  Axisymmetric  Boundary  Conditions 

For  the  given  flow  fleld,  there  are  flve  boundaries,  where  boundary  conditions  have  to  be  specified. 

Inflow 

At  the  inflow  boundary  all  flow  variables,  except  for  the  pressure  within  the  subsonic  layer, 
are  specified,  given  by  a  flat  plate  boundary  layer  similarity  solution.  In  the  subsonic  layer 
the  pressure  at  the  inflow  boundary  is  assumed  to  be  constant  in  streamwise  direction  and 
is  extrapolated  from  the  flow  field,  using 


Pl,n=P2,n,  if  Afi,n  <  1- 


(3.28a) 


Outflow 

At  the  outflow  boundary  the  second  derivative  in  the  axial  direction,  is  set  to  zero. 
This  leads  to  the  following  fourth-order  accurate  one-sided  finite  difference  approximations 
for  the  first  derivatives  near  the  outflow  boundary: 


3A(f>M,n  +  “  42<^M-2,n  +  50M-3,n 

(3.28b) 

M—l,n 

66Ae 

d(f) 

8^4>M,n  —  1O80M-I,n  +  T7(f)M-2,n  — 

(3.28c) 

M,n 

66Ae 

Walls 

At  the  walls,  all  velocities  are  set  to  zero.  For  the  calculation  of  a  steady  flow  field,  the 
body  is  assumed  to  be  adiabatic.  Therefore  the  normal  gradient  of  the  temperature  is 
zero.  For  the  calculation  of  an  unsteady  flow  field,  the  wall  temperature  is  kept  constant 
in  time.  The  pressure  at  the  wall  is  extrapolated  from  the  flow  field,  using  the  wall  normal 
momentum  equation.  Thus,  the  conditions  are 
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At  the  horizontal  wall: 

'^r\m,w  — 

Tm.w  =  [300  (Tm, w+1  —  Tm,w+2)  +  200Tnj,u;+3  —  75Tm,w+4  +  12T,n,w+5]  > 

for  a  steady  flow  fleld, 

Tm,w  =  ,  for  an  unsteady  flow  field, 

Pm,w  =  ~Trrm,w  +  [300  {pm,w+l  +  Pm,'w+'i.'^'m,w+\  "I"  ’^rrm,w+\  ~  Pm, w+2 

~  Pm,w+2'^m,w+2  ~  ‘^rrm,w+^  d"  (pm,w+2  d"  Pm,w+3f^m,w+3 
d'T'rrmjM+a)  ~  75  (pTn,w+4  d-  Pm,‘w+4v‘m,'w+4  d"  Trrm,w+^ 
d"12  (^Pm,w+5  d"  Pm,w+5^m,w+5  d*  TrrmjW+s)] 

^  i'^rxm-2,w  “ 


•2,11)  ^'^rxm—l,w  d"  STVim+l,!!)  '^rxm+2,w) 


At  the  base: 

^*l6,n  ~  0) 

^r\t,n  =  0) 

Tb,n  =  [300  {Tb+I^n  —  76+2, n)  d-  200T6+3,n  —  75r6+4_„  +  12Tb+5^„]  , 

for  a  steady  flow  field, 

Tb,n  =  Tolj,  „ ,  for  an  unsteady  flow  field, 

Pb,n  =  ~Txxb,n  d"  [300  (p6+l,n  d-  Txxb+l,n  ~  Pb+2,n  ~  Txxb+2,n) 

+200  (p6+3,n  d"  Txxb+3,n)  75  {pb+4,n  d"  '^ii6+4,n) 

+  12  {pb+5,n  d"  7116+5,11)]  ~  60A^  — 

d-Srrift  „+i  -  Trxb,n+2)  d"  ^  •  (3.28e) 

np)  J  6 

Axis 

At  the  axis  the  symmetry  of  the  flow  field  is  used  to  determine  the  flow  quantities  at 
locations  outside  the  integration  domain.  Therefore,  it  follows 

P\m,--\  ^lm,2  i  P\m,—2  P\m,3  > 

P^x\m,-l  —  P^®l77i,2  5  P^a:lm,-2  P^x\m,3  ) 


L  12An  ^-^rxb, n-l 
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P^r\m,-X  =  -  P^rU,2  I  PVr\m,-2  =  “  P^r\m,Z  . 

P^\m-l  =  Pelm.2  ;  P^\m,-2  =  P^\m,Z  ■  (3-28f) 

For  all  ^  terms  the  rule  of  THopital  is  applied  as  described  in  appendix  C. 

Free  stream 

At  the  free  stream  boundary,  three  different  boundary  conditions  are  applied,  according 
to  the  free  stream  Mach  number  of  the  flow  field.  The  conditions  axe  described  in  detail 
in  Appendix  B. 

Initial  condition 

As  initial  condition  for  the  calculation  of  a  steady  axisymmetric  flow  field  the  similarity 
solution  of  a  compressible  flat  plate  boundary  layer  is  used.  In  addition,  the  radial  velocity 
is  determined  by  integration  of  the  steady  continuity  equation.  For  the  calculation  of  an 
unsteady  axisymmetric  flow  field,  a  previously  calculated  steady  axisymmetric  flow  field 
is  used  as  initial  condition. 

3.5.2  Three  Dimensional  Boundary  Conditions 
Inflow 

For  the  three-dimensional  (unsteady)  calculations,  all  variables  at  the  inflow  boundary 
are  specified.  The  only  exception  is  the  pressure  in  the  subsonic  layer,  where  the  first 
derivative  of  the  pressme  is  set  to  zero.  So, 

Pi,n  =  P2,n  if  <  1  .  (3.29a) 


Outflow 

As  for  the  axisymmetric  calculations,  the  second  derivatives  of  all  flow  quantities  are  set 
to  zero  at  the  outflow  boundary.  For  the  first  derivatives,  the  same  difference  formulas  as 
for  the  axisymmetric  flow  are  used. 

Walls 

At  the  walls,  all  velocities  are  set  to  zero.  The  temperature  along  the  walls  is  kept  constant 
in  time,  assuming  the  body  cannot  change  temperature  on  the  given  time  scale.  The  wall 
pressure  is  extrapolated  in  physical  space  in  the  same  fashion  as  for  the  axisymmetric 
problem.  This  leads  to  the  following  boundary  conditions  at  the  azimuthal  location  0  = 
hAO: 

At  the  horizontal  wall: 

'^r\m,w,h  ~ 

'^s\m,w,h 
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Tm,w,h  — 

Pm,w,h  —  ~'^rrm,w,h  ^300  {^rn,w+X,h  +  Pm,w+\,hVm,w-Xl,h  '^rrm,tu+l,/i 

~Pm,iu+2,/i  ■“  Pm,,w-i-2,h'^m,w+2,h  ~  '^rrm,u;+2,/i) 

+200  (prn,w+3,h  +  Pm,«)+3,ft^m,w+3,/i  '^rrm,w+Z,l^ 

—  75  (pm,w+A,h  +  PTn,w+4,hf^m,w+A,h  '^rrm,tt)+4,/i) 

+12  (pm.,w+5,h  +  Pm,w+5,h‘f^m,w+5,h  "I"  '^rrm,tu+5,/i)] 


(3.29b) 


dr 

dr] 

W 

,dx 

{'^rxm—2,v},h  ^'^rxm—l,w,h 


,m  12Ae 

A8Tj-XTn+l,w,h  '^rxm+2,w,h)  "b  '^rrm,w,h  '^06m,w,h  "b 


dTffr 


de 


m,w,h 


m,w,h 


is  calculated  in  Fourier  space  and  then  trans- 


where  the  azimuthal  derivative  ^ 
formed  into  physical  space. 

At  the  base: 


'^r\b,n,h  — 

Mb,n,h  =  0i 
Tb,n,h  =  Tolj  , 

Pb,n,h  —  ~'^xxb,n,h  +  [300  (p6+i,n,ft  +  Txxb+l,n,h  ~  Pb+2,n,h  ~  Txxb+2,n,h) 


+200  {pb+3,n,h  +  Txxb+3,n,h)  75  {pb+4,n,h  +  '^xxb+4,n,h.) 


+  12  {pb+5,n,h  +  Txxb+b,n,h)]  “  60A^ 


dx 


dr] 

b 

.  dr 

nl2Ar? 


ijrxb,n— 


2,h 


^'^rxb,n—lh  +  ^'^rxb,n+X,h  '^TXb,n+2,h) 


.  1  1  _  .  9'^ex 

r{v)  V  d9 


b,n,h  ) 


(3.29c) 


Axis 

The  Fourier  expansion  of  the  three-dimensional  governing  equations  results  in  an  addi¬ 
tional  boundary  at  the  axis  of  symmetry.  Therefore,  a  set  of  boundary  conditions  has 
to  be  described  for  the  discretized  equations.  At  the  axis,  the  even  Fourier  modes  are  of 
symmetric  nature  and  the  odd  modes  are  of  antisymmetric  nature.  They  can  be  prescribed 
at  the  boundary  as  follows.  The  derivation  of  the  boundary  conditions  for  the  different 
Fourier  modes  is  described  in  detail  in  appendix  C. 

fc  =  0 

The  boundary  conditions  for  Fourier  mode  k  =  0  are  exactly  the  same  as  for  the 
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axisymmetric  case, 


^Im,— 1,0  ^lm,2,0 

>  ^Im,— 2,0  ^lm,3,0  ’ 

!  ^3:|,n,-2,0  =  P^x\m,'i,0-: 

“  P^r\m, 2,0 

i  Ini, -2,0  =  ~  P^r|m,3,0  > 

P^\m,—  \,0  P^\m, 2,0 

)  P^\m, -2,0  ~  ^^lm,3,0  > 

(3.29d) 

while  for  fc  =  1 


=0> 

=  0, 

and  for  /c  >  2 


P^m,X,k 

=  0, 

P^x\m,\,k 

=  0, 

P^^\m,\,k 

=  0, 

=  0, 

P^\m,\,k 

=  0. 

For  all  ^  terms  the  rule  of  I’Hopital  is  applied  as  described  in  appendix  C. 


(3.29e) 


(3.29f) 


Free  stream 

As  for  the  axisymmetric  flow,  three  difierent  types  of  boundary  conditions  are  applied  at 
the  free  stream  boundary,  according  to  the  free  stream  Mach  number.  These  are  described 
in  more  detail  in  appendix  B. 


Initial  conition 

As  initial  condition  a  previously  calculated  steady  axisymmetric  flow  fleld  is  used. 


3.6  Filtering 

For  the  Large-Eddy  Simulations  the  flow  quaintities  p,  p^Vi  and  We  are  Altered  at  every 
Runge-Kutta  stage.  In  the  present  research  a  sixth-order  compact  filter  has  been  chosen  [see, 
for  example,  [Lele  (1992)]].  This  Alter  has  been  tested  extensively  by  [Bachman  (1996)].  As  an 
example,  the  equation  for  a  quantity  filtered  in  the  axial  direction  at  a  given  location  (xi)  is 

PFfi-2  +  Olpfi-l  +  fi  +  “f/x+I  +  PFfi+2 

=  o/i  +  2  (-^*-1  +  2  (-^*-2  +  /i+2)  +  2  -^*+3)  (3.30) 

where  ft  are  the  unfiltered  quantities  and  /j  are  the  filtered  quantities. 
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The  coefficients  are  given  by 

ap  =  0.6522474  /3f  =  0.1702929 

a  =  (2  +  3aF)/4  6  =  (9  +  160^  +  10/5f)/16  (3.31) 

c={aF  +  4M/^  rf=(6/3F-l)/16 

This  filter  has  been  successfully  used  for  Large-Eddy  Simulations  of  incompressible  fiat  plate 
boundary  layers  by  [Bachman  (1996)]. 

At  all  Reynolds  numbers  calculated  in  the  present  research,  the  boundary  layer  approaching 
the  corner  of  the  base  remained  laminar  and  turbulence  modeling  was  not  needed  in  that  region. 
The  filtering,  thus,  is  applied  solely  in  the  second  subdomain,  downstream  of  the  base  of  the 
axisymmetric  body.  In  the  transition  region  between  the  two  domains,  the  turbulence  model  is 
slowly  ramped  in,  using  the  same  formula  as  for  the  wall  damping  function  (see  section  2.4.7). 

3.7  Disturbance  Generation 

For  the  unsteady  (distmrbed)  calculations,  disturbances  are  introduced  into  the  flow  field.  In 
general,  there  are  two  difierent  types  of  disturbances  that  can  be  generated,  pulse  disturbances 
and  continuous  disturbances.  The  main  difierence  between  the  two  types  of  disturbances  lies 
in  their  frequency  spectra.  A  continuous  disturbance  typically  has  a  fixed  frequency  and  thus, 
only  an  isolated  single  frequency  disturbance  is  introduced  into  the  flow  field  (see  figure  3.3). 
This  allows  the  investigation  of  the  stability  behavior  for  one  given  frequency. 


1.0 

0.5 

0.0 

0.0  0.2  0.4  0.6  0.8  1.0 

frequency 

Figure  3.3:  Continuous  disturbance;  time  signal  (left)  and  Fourier  transform  (right). 

On  the  other  hand,  a  pulse  disturbance  generates  a  broad  spectrum  of  frequencies  in  the 
flow  field.  The  broadness  of  the  spectrum  depends  on  the  length  of  the  pulse  (see  figure  3.4 
and  figure  3.5).  Thus,  the  broadest  spectrum,  or  in  other  words  the  widest  frequency  band,  is 
generated  by  a  single  spike  disturbance  in  time  (see  figure  3.6). 
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time  frequency 


Figure  3.4:  Single  pulse  disturbance  over  three  periods;  time  signal  (left)  and  Fourier  transform 
(right). 


Figure  3.5:  Single  pulse  disturbance  over  one  period;  time  signal  (left)  and  Fourier  transform 
(right). 

3.7.1  Continuous  Disturbances 

Disturbances  of  different  Fourier  modes  are  introduced  into  the  flow  fleld  by  radial  blowing  and 
suction  through  a  circular  slot  near  the  base  of  the  axisymmetric  body  (see  figure  3.7).  When  a 
normal  velocity  distribution  of  the  form  shown  in  figure  3.7  is  used,  the  net  mass  flow  through  the 
disturbance  slot  is  zero  at  every  instant  of  time.  This  technique  produces  predominantly  vorticity 
disturbances  and  only  to  a  lesser  extent  undesirable  ’’sound”  disturbances.  The  disturbances 
first  develop  within  the  boundary  layer  region  and  then  travel  into  the  free  shear  layer  region 
and  the  recirculation  zone. 

By  varying  the  real  and  imaginary  parts  of  the  Fourier  component  of  the  disturbed  radial 
velocity  at  the  blowing  and  suction  slot,  rotating  disturbances  can  be  introduced,  thus  pro¬ 
ducing  a  helical  wave  that  travels  downstream.  It  has  been  shown  for  incompressible  wake 
flows  that  the  flow  field  is  unstable  with  respect  to  these  helical  disturbances  [see,  for  example, 
[Schwarz  (1996)]]. 
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time  frequency 

Figure  3.6:  Spike  disturbance;  time  signal  (left)  and  Fourier  transform  (right). 

3.7.2  Pulse  Disturbances 

A  pulse  disturbance  can  also  be  introduced  through  the  blowing  and  suction  slot,  by  varying  the 
amplitude  over  time.  This  generates  a  wave  packet  and  allows  investigation  of  the  disturbance 
behavior  for  a  spectrum  of  frequencies  as  the  wave  packet  travels  downstream. 

In  addition  to  blowing  and  suction,  a  second  method  is  used  for  the  generation  of  pulse 
disturbances.  In  this  case,  a  discrete  pulse  is  introduced  locally  into  the  flow  fleld  at  a  location 
within  the  recirculation  region.  Usually,  a  pulse  in  form  of  a  single  spike  in  the  flrst  azimuthal 
Fourier  mode  of  the  density  at  a  specifled  grid  point  and  only  one  instance  in  time  is  used.  This 
method  generates  disturbances  with  the  broadest  possible  spectrum  of  frequencies.  It  is  mainly 
used  to  investigate  the  existence  and  behavior  of  absolute  instabilities. 

3.8  Computational  Procedure 

For  the  investigation  of  the  behavior  of  three  dimensional  flow  disturbances,  an  axisymmetric 
initial  undisturbed  flow  field  has  to  be  determined.  Since  no  similarity  solution  exists  for  the 
compressible  axisymmetric  wake,  the  first  step  in  the  computational  procedure  is  to  calculate  a 
steady  flow  field  using  an  unsteady  axisynunetric  Navier-Stokes  program.  As  initial  condition 
for  this  axisymmetric  calculation,  a  similarity  solution  for  a  compressible  flat  plate  boundary 
layer  is  used.  The  imsteady  calculation  is  terminated  once  the  temporal  changes  of  the  flow 
variables  are  much  smaller  than  the  expected  disturbance  amplitudes. 

This  axisymmetric  flow  field  is  then  used  as  initial  condition  for  the  three  dimensional  dis¬ 
turbance  calculation  .  Disturbances  are  introduced  either  through  a  blowing  and  suction  slot 
along  the  circular  body  near  the  base  or  locally  into  the  flow  field  as  a  single  pulse  of  width  At 
(see  section  3.7).  The  time  dependent  evolution  of  the  disturbances  is  monitored  over  several 
time  periods  to  investigate  growth  or  decay  and  influences  on  the  global  flow  field. 
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Figure  3.7:  Disturbance  generation  through  a  blowing  and  suction  slot. 

For  suflBciently  high  Reynolds  numbers,  local  regions  of  very  high  velocity  gradients  appear 
in  the  flow  fleld.  In  that  case,  the  calculation  with  DNS  is  terminated  before  the  gradients 
become  so  large  that  they  would  cause  numerical  instability  and  eventually  a  termination  of 
the  program.  The  instantaneous  flow  field  at  this  point  is  used  as  an  initial  condition  for  a 
Large-Eddy  Simulation. 


Chapter  4 


Code  Validation 


At  the  present  time,  no  data  are  available  for  supersonic  axisymmetric  waJce  flow  flelds  in  the 
Reynolds  number  range  between  Brd  —  500  and  ifep  =  30, 000.  Therefore,  most  of  the  code  val¬ 
idation  was  performed  by  comparing  the  results  of  very  low  Mach  number  subsonic  calculations 
to  DNS  data  for  an  incompressible  axisymmetric  wake  that  was  obtained  by  [Schwarz  (1996)] 
[see  also  [Schwarz  et  al.  (1994)]]  and  by  further  comparison  to  experimental  results  that  were 
obtained  in  our  own  water  channel  by  [Siegel  (1996)].  In  addition,  for  supersonic  flow  a  step 
size  analysis  was  performed  for  the  axisymmetric  flow  case. 

4.1  Axisymmetric  Flow 

As  a  first  step,  the  code  validation  was  performed  for  steady  and  unsteady  axisymmetric  flow 
calculations.  The  code  validation  was  done  for  a  low  Mach  number  subsonic  flow  by  comparison 
with  results  from  incompressible  calculations  and  from  experiments.  For  supersonic  flow  a  step 
size  analysis  was  performed,  because  of  the  lack  of  experimental  data  and  previous  numerical 
results. 


4.1.1  Subsonic  Flow  Fields 

For  validation  of  the  Navier-Stokes  code  for  the  steady  (undisturbed)  axisymmetric  flow,  the  flow 
field  for  a  free  stream  Mach  number  of  Moo  =  0.2  and  a  global  Reynolds  number  oi  Ben  =  Ij  000 
was  calculated  for  two  difierent  spatial  resolutions.  In  both  cases  the  grid  was  stretched  in 
the  axial  and  radial  directions  with  the  same  stretching  factor.  The  results  were  compared  to 
results  from  a  calculation  for  an  incompressible  axisymmetric  waJce  by  [Schwarz  (1996)].  For  the 
incompressible  calculations  no  grid  stretching  was  used.  The  parameters  for  the  incompressible 
calculation  by  Schwarz  (1996)  and  cases  Al  and  A2  of  the  current  compressible  calculations  are 
listed  in  table  A.l. 

For  comparisons  of  the  compressible  and  the  incompressible  simulations,  the  length  of  the 
recirculation  zone  was  chosen  as  a  global  characteristic  of  the  steady  flow  field,  and  the  azimuthal 
vorticity  at  the  corner  of  the  base  was  chosen  as  a  local  characteristic.  The  recirculation  length. 
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which  is  the  distance  of  the  free  stagnation  point  from  the  base,  for  the  current  calculation  is 

Lrec  =  2.1,  (4.1) 

for  all  three  cases.  Thus,  both  compressible  calculations  compare  very  well  with  the  result  of  the 
incompressible  calculations.  Figures  4.1  and  4.2  show  the  steady  streamlines  for  both  validation 
cases.  The  streamlines  clearly  exhibit  the  region  of  recirculation.  Prom  the  corner  of  the  base, 
the  flow  separates  geometrically  and  bends  gradually  towards  the  axis  of  symmetry.  The  free 
stagnation  point  can  be  found  at  re  =  2.1  for  both  grid  resolutions.  Both  cases  show  basically 
identical  streamlines  and  are  in  very  good  agreement  with  the  incompressible  results. 


n  I  I  \  i  I - 1 - 1 - 1 — — - 1  i  r- 
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Figure  4.1:  Steady  streamlines  for  validation  case  Al. 
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Figure  4.2:  Steady  streamlines  for  validation  case  A2. 
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The  azimuthal  vorticity  at  the  corner  of  the  base,  which  is  defined  as 

(jj  =  —V  X  u,  (4.2) 


turns  out  to  be  dependent  on  the  grid  resolution  at  the  corner  in  both  the  incompressible  and  the 
compressible  calculations.  At  the  corner  one-sided  differences  have  been  used  for  approximation 
of  the  spatial  derivatives  of  the  velocity.  For  the  incompressible  calculation,  corner  vorticity  was 
found  to  be 


«.incomp|^^o  ,.^0  5 


=  6.8, 


whereas  for  cases  A1  and  A2  the  values  are 


‘^fl.iL=o,r=o.5  =  7-2,  and  a;e,2L=o,r=o.5  =  ^.6.  (4.4) 

Considering  the  different  resolution  in  the  three  cases  and  the  different  methods,  the  values  are 
in  reasonable  agreement.  Figures  4.3  and  4.4  show  isolines  of  azimuthal  vorticity  for  the  two 
compressible  validation  cases.  In  both  plots,  the  vorticity  levels  for  the  isolines  are  identical 
duel  =  1,2,3...).  The  lines  in  the  two  cases  are  almost  identical  and  also  coincide  with  the  plot 
shown  in  Schwarz  (1996).  In  conclusion,  the  agreement  for  the  steady  axisymmetric  fiow  for  the 
case  of  a  very  low  subsonic  free  stream  Mach  number  is  very  good. 


X 

Figinre  4.3:  Steady  azimuthal  vorticity  field  for  validation  case  Al. 
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For  validation  of  the  Navier-Stokes  code  for  the  unsteady  (disturbed)  axisymmetric  flow, 
the  response  of  the  subsonic  flow  fleld  of  case  A1  (obtained  from  the  steady  calculation  as  dis¬ 
cussed  above)  to  a  small  axisymmetric  distm-bance  was  calculated  (case  A3).  In  this  calculation, 
the  steady  axisymmetric  flow  was  disturbed  with  a  continuous  axisymmetric  disturbance  at  a 
frequency  of  /3  =  1.0.  The  disturbance  was  generated  through  a  blowing  and  suction  slot  as  de¬ 
scribed  in  section  3.6.  The  amplitude  of  the  blowing  and  suction  distribution  was  0.001  percent 
of  the  free  stream  velocity.  Therefore,  a  linear  behavior  of  the  generated  disturbance  wave  can 
be  expected. 

The  calculation  was  carried  out  until  the  entire  flow  field  reached  a  time  periodic  state.  At 
this  point  the  spatial  development  of  the  disturbance  waves  in  the  wake  region  was  investigated 
to  obtain  the  spatial  growth  rates.  Figure  4.6  shows  the  spatial  growth  rate  of  the  distmbance 
frequency  versus  downstream  location.  For  comparison  with  the  incompressible  calculations  the 
disturbance  waves  were  monitored  in  the  form  of  the  square  root  of  the  kinetic  disturbance 
energy.  Assuming  a  behavior  of  the  disturbance  waves  of  the  form 

V^  =  (4.5) 

with  complex  a  and  /3,  the  spatial  growth  rate  is  given  by  the  imaginary  part  of  the  wavenumber 
as  follows 


The  results  show  a  region  of  positive  growth  (denoted  by  a  negative  Oj),  which  reaches  up 
to  about  1.2  diameters  downstream  of  the  base.  Downstream  of  this  location  the  disturbances 
are  damped  and  eventually  die  out.  Comparing  these  growth  rates  to  linear  stability  theory 
and  incompressible  calculations  [Figure  4.4  in  Schwarz  (1996)]  very  good  agreement  was  found 
(see  figure  4.5).  Thus,  also  for  the  unsteady  axisymmetric  calculations  very  good  agreement  was 
found  between  the  compressible  and  the  incompressible  calculations. 

4.1.2  Supersonic  Flow  Fields 

As  mentioned  earlier,  no  experimental  data  are  available  for  supersonic  wakes  in  the  Reynolds 
number  regime  that  is  applicable  for  Direct  Numerical  Simulations.  For  that  reason  a  step- 
size  investigation  was  performed  for  an  axisymmetric  flow  at  a  free  stream  Mach  number  of 
Moo  =  2.46  and  a  global  Reynolds  number  of  Ebd  =  30,000.  The  flow  field  was  calculated  for 
two  different  spatial  resolutions.  The  computational  parameters  for  both  cases  (A4  and  A5)  are 
given  in  table  A.l. 

Comparisons  were  made  for  the  axial  and  radial  velocities,  the  density,  the  pressure,  the 
temperature  and  the  local  Mach  number.  The  results  in  form  of  isolines  are  shown  in  figures 
4.6  through  4.11,  respectively.  In  all  the  figures,  the  results  for  the  lower  resolution  (case  A4) 
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Figure  4.5:  Spatial  growth  rate  ai  for  an  axisymmetric  disturbance  of  =  1.0,  /3j  =  0.0; 
Red  =  1,000  [circles  from  Linear  Stability  Theory  by  Schwarz  (1996)]. 


are  shown  at  the  top  and  for  the  higher  resolution  (case  A5)  below.  All  flow  quantities  show 
excellent  agreement  between  the  two  different  cases.  The  wiggles  at  some  isolines  for  the  pressure 
in  the  lower  resolved  case  (see  figures  4.7  through  4.9)  are  due  to  the  plotting  software  (linear 
interpolation) .  These  results  indicate  that  the  lower  resolution  was  sufiicient  for  the  calculation 
of  this  flow  field.  Coincidentally,  this  resolution  was  the  lowest  possible  to  avoid  numerical 
instabilities  in  the  calculations. 


Figure  4.6:  Isolines  of  axial  velocity  at  Moo  ^  2.46  and  Ren  =  30,000  for  two  different  spatial 
resolutions  [Cases  A4  (top)  and  A5  (below)]. 
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Figure  4.8;  Isolines  of  density  at  Mqo  =  2.46  and  ife/j  =  30,000  for  two  different  spatial 
resolutions  [Cases  A4  (top)  and  A5  (below)]. 
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resolutions  [Cases  A4  (top)  and  A5  (below)]. 


X 

Figure  4.10:  Isolines  of  temperature  at  Moo  =  2.46  and  Rbd  =  30,000  for  two  different  spatial 
resolutions  [Cases  A4  (top)  and  A5  (below)]. 
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4.2  Three-dimensional  Flow 

A  more  thorough  validation  would  be  desirable  at  this  point  for  the  codes  used  for  the  calcu¬ 
lation  of  three-dimensional,  time  dependent  supersonic  flow.  However,  the  three-dimensional 
simulations  at  the  chosen  resolution  for  the  supersonic  flows  were  already  very  CPU  time  and 
memory  intensive  (see  chapter  5).  A  three-dimensional  calculation  at  the  given  resolution  of 
case  A4  typically  required  600  MBytes  of  main  memory  and  in  the  order  of  1,000  CPU  hours 
on  a  CRAY  C90.  Therefore,  the  three-dimensional  code  was  validated  for  a  subsonic  flow  at 
Moo  =  0.2  and  Reo  =  1,000. 

For  the  validation  of  the  Navier-Stokes  code  for  calculating  the  three-dimensional  unsteady 
(disturbed)  flow,  the  response  of  the  subsonic  flow  field  of  case  A1  (obtained  from  the  steady 
calculation  as  discussed  above)  to  a  three  dimensional  disturbance  input  was  calculated  and  the 
results  were  compared  to  results  from  incompressible  calculations  by  [Schwarz  (1996)]  and  to 
water  channel  experiments  carried  out  in  the  Hydrodynamic  Laboratory  at  the  University  of 
Arizona.  The  simulation  was  performed  with  up  to  eight  complex  Fourier  modes  in  azimuthal 
direction  and  the  same  spatial  resolution  as  for  case  Al. 

4.2.1  Absolute  Instability 

For  a  first  calculation  (case  Tl),  the  flow  was  disturbed  with  a  pulse  disturbance  of  the  frequency 
/?  =  0.1  through  the  blowing  and  suction  slot.  The  pulse  had  a  duration  of  two  full  periods  and 
was  ramped  up  and  down  in  time,  using  a  sin^  function.  In  a  second  calculation,  the  flow  was 
disturbed  only  at  the  very  first  time  step  locally  within  the  recirculation  region  (as  described  in 
section  3.6)  at  2C  =  0.25  and  r  =  0.25.  The  flow  response  is  similar  for  both  the  compressible 
and  the  incompressible  calculations. 

In  incompressible  simulations,  [Schwarz  (1996)]  found  that  for  the  global  Reynolds  number 
ifejo  =  1, 000  the  flow  field  includes  a  region  of  absolute  instability  with  regard  to  helical  distur¬ 
bances  of  the  first  azimuthal  Fourier  mode.  Figure  4.12  shows  the  time  response  of  the  radial 
momentum  at  the  centerline  in  the  near  waJce,  half  a  radius  downstream  of  the  base.  In  this 
case,  the  disturbance  was  generated  through  blowing  and  suction.  The  response  clearly  shows  an 
exponential  growth  of  the  distmbance,  which  indicates  the  existence  of  an  absolute  instability, 
as  has  been  found  in  the  incompressible  simulations. 

For  a  second  calculation  (case  T2),  a  disturbance  was  generated  by  a  single  pulse  in  the 
near  wake  region,  as  explained  before.  Figure  4.13  shows  the  time  response  of  the  flow  in  the 
form  of  the  radial  momentum  of  four  difierent  azimuthal  Fourier  modes  at  the  location  of  the 
disturbance  input.  In  this  case  the  flow  starts  to  deviate  from  the  axisymmetric  state  and 
the  signal  keeps  oscillating  in  time  in  a  non-periodic  fashion.  This  was  also  observed  for  the 
incompressible  calculations.  From  this  point,  quantitative  comparison  is  not  very  meaningful, 
since  the  fluctuations  are  not  periodic. 

Therefore,  the  comparison  of  the  instantaneous  flow  field  with  experiments  and  incompress¬ 
ible  simulations  can  only  be  of  a  qualitative  nature.  For  this  qualitative  comparison,  particles 


CHAPTER  4.  CODE  VALIDATION 


55 


Figure  4.12:  Time  response  of  radial  momentum  at  r=0  to  a  three  dimensional  blowing  and 
suction  pulse  disturbance  in  the  near  waJce  at  Mqo  =  0.2  and  ifep  =  1,000. 

have  been  released  into  the  flow.  The  method  is  described  in  section  5.1.1.  Figure  4.14  shows 
side  views  of  the  instantaneous  flow  fleld.  The  views  are  perpendicular  to  each  other.  They  show 
that  the  flow  is  symmetric  to  the  x-y  plane  (0  =  0).  This  phenomenon  could  also  be  observed 
in  experiments  and  incompressible  simulations  as  presented  by  Schwarz. 

Comparing  the  structures  to  pictmres  shown  in  Schwarz  (1996)  a  good  similarity  in  the  shape 
of  the  structures  is  seen.  For  further  simulations  this  symmetry  was  prescribed,  so  that  only 
half  of  the  points  in  physical  space  need  to  be  calculated  and  main  memory  and  CPU  time  are 
saved.  In  figure  4.15,  cuts  through  the  y-z  plane  at  four  different  x  locations  are  presented. 
The  plots  show  structures  very  similar  to  the  structures  that  have  been  observed  in  water 
channel  experiments  by  [Siegel  (1994)]  [see,  for  example,  figure  4.16].  To  further  investigate  the 
development  of  the  disturbances  and  allow  a  better  comparison  with  the  experimental  results, 
the  flow  was  disturbed  by  a  periodic  disturbance  [see  section  4.2.2]. 
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Figure  4.13:  Time  response  of  radial  momentum  to  a  single  pulse  of  Fourier  mode  A:  =  1  at 
r=0.25  and  z=  0.25  at  Mqo  =  0.2  and  Red  =  1,000. 
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Figure  4.14;  Flow  visualization  via  paxticles  at  Moo  =  0.2  and  Reo  =  1,000,  viewed  from 
0  =  7r/2  (top)  and  from  0  =  0  (bottom). 
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Figure  4.15:  Flow  visualization  via  particles  at  Moo  =  0.2  and  ifep  =  1,000,  cross-sectional  cut 
at  I  =  2.0. 


Figure  4.16:  Flow  visualization  via  particles  for  ife/j  =  1,000,  cross-sectional  cut  at  a:  —  2.0 
(water  channel  experiment  by  Siegel  (1994). 
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4.2.2  Periodic  Disturbances 

In  order  perform  a  quantitative  comparison  with  the  experimental  results  the  same  flow  fleld 
as  above  (case  Al)  was  disturbed  continuously  with  a  flxed  disturbance  frequency  (case  A3). 
For  the  generation  of  the  periodic  disturbances,  the  blowing  and  suction  slot  was  used.  The 
disturbances  consisted  purely  of  azimuthal  modes  A:  =  —  1  and  k  =  1,  resembling  two  counter¬ 
rotating  disturbances  of  the  first  helical  Fourier  mode.  The  Strouhal  number  of  the  disturbance 
based  on  the  free  stream  velocity  and  the  base  diameter  was  chosen  to 

StD  =  (^  =  0.159.  (4.8) 

Uqo 

A  flrst  attempt,  using  a  disturbance  amplitude  of  Ai^i  =  0.01  did  not  show  any  significant 
influence  of  the  disturbance  on  the  flow  fleld.  Therefore,  in  a  second  and  flnal  calculation  (T3), 
a  disturbance  amplitude  of 

Ai,i  =  0.1  (4.9) 

was  used.  The  amplitude  resembles  the  maximum  blowing  and  suction  velocity,  based  on  the 
free  stream  velocity. 


Figure  4.17:  Time  response  of  total  vorticity  to  a  periodic  disturbance  of  Fourier  mode  fe  =  1 
with  StD  =  0.159  and  A  =  0.1  at  Moo  =  0.2  and  Bbd  =  1,000  {x  =  0.25,  r  =  0.25). 

The  local  time  response  of  the  total  vorticity  at  a  point  within  the  recirculating  region 
{x  =  0.25,  r  =  0.25)  is  shown  in  figure  4.17.  The  generation  of  a  small  amplitude  disturbance 
(Ai,i  =  0.01)  was  started  at  the  beginning  of  the  calculation  {t  =  70).  As  mentioned  earlier,  the 
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flow  did  not  show  any  visible  reaction  to  the  disturbance  (see  figure  4.17).  The  signal  remains 
nori-periodic.  At  a  later  stage  {t  =  250)  the  amplitude  was  raised  to  the  final  value.  After 
several  time  periods  the  response  of  the  flow  becomes  periodic  with  the  same  frequency  as  the 
generated  disturbance  (see  figure  4.17).  At  this  point,  the  structures  also  appear  in  a  periodic 
fashion. 

For  a  first  validation,  the  periodic  flow  was  calculated  over  two  time  periods  for  three  different 
spatial  resolutions.  The  different  parameters  are  shown  in  table  A.2  for  the  cases  T3  through 
T5.  The  results  for  the  three  different  cases  were  Fourier  transformed  in  time  to  separate  the 
fundamental  disturbance  frequency  and  its  higher  harmonics.  Figures  4.18  through  4.21  show 
the  resulting  velocity  fields  for  the  several  Fourier  modes  in  time  and  azimuthal  direction  for 
the  three  different  cases.  It  appears  that  they  all  compare  very  well  in  the  region  immediately 
downstream  of  the  base.  Further  downstream  the  results  for  the  different  cases  show  small  but 
noticeable  differences.  This  is  due  to  the  grid  stretching,  which  results  in  a  lower  resolution  in 
the  region  further  downstream.  This  inaccuracy,  however,  does  not  seem  to  affect  the  solution 
in  the  near  wake  region.  Therefore,  for  further  comparison  with  the  experiments  the  lower 
resolution  case  was  used. 

Because  of  the  lower  resolution,  all  simulations  could  be  performed  on  a  single  processor  of 
a  Silicon  Graphics  Power  Challenge  at  the  Computational  Fluid  Dynamics  Laboratory  of  the 
Aerospace  and  Engineering  Department  at  the  University  of  Arizona.  A  typical  simulation  used 
about  20  CPU  hours  for  the  calculation  of  one  time  period. 
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Figure  4.18:  Isolines  of  the  mean  of  axial  -selocity  for  three  different  spatial  resolutions  for 


Moo  =  0.2  and  Beo  =  1,000  [cases  T3  (top)  through  T5  (bottom)]. 


CHAPTER  4.  CODE  VALIDATION 


61 


Figure  4.19:  Isolines  of  the  fundamental  of  first  azimuthal  mode  of  axial  velocity  for  three  dif¬ 
ferent  spatial  resolutions  for  Mqo  =  0.2  and  ife/j  =  1,000  [cases  T3  (top)  through  T5  (bottom)]. 

In  the  experiment,  which  was  conducted  in  the  water  channel  at  the  Hydrodynamics  Lab¬ 
oratory  of  the  Aerospace  and  Mechanical  Engineering  Department  by  [Siegel  (1996)],  the  fiow 
was  disturbed  using  a  vibrating  ribbon.  The  Strouhal  number  was  the  same  as  for  the  numerical 
simulation.  However,  the  amplitude  of  the  disturbance  in  the  experiment  was  not  determined. 
As  in  the  numerical  simulation,  the  flow  becomes  periodic  after  several  periods  of  the  generated 
disturbance.  In  addition,  the  structmes  appeared  to  remain  symmetric  to  the  symmetry  plane 
of  the  vibrating  ribbon.  Thus,  the  plane  of  symmetry  was  fixed,  in  contrast  to  the  self  excited 
axisymmetric  wake. 

For  a  comparison  between  the  simulations  and  the  experiments,  velocity  measurements  were 
made  at  a  cross  section  which  was  about  3.5  diameters  downstream  of  the  blunt  base.  For 
the  measTuements  a  hot  film  anemometer  was  used.  This  allowed  only  a  measurement  of  the 
combination  of  two  velocity  components.  The  hot  film  probe  was  directed  such  that  in  the  plane 
of  symmetry  it  would  measure  the  combination  of  axial  and  radial  velocity.  Thus,  results  were 
compared  for  the  combined  quantity 

q  =  (4.10) 

The  measiued  data,  as  well  as  the  calculated  results,  were  Fourier  transformed  in  time  in 
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Figure  4.20:  Isolines  of  the  fundamental  of  first  azimuthal  mode  of  radial  velocity  for  three  dif¬ 
ferent  spatial  resolutions  for  Moo  =  0.2  and  Beo  =  1,000  [cases  T3  (top)  through  T5  (bottom)]. 

order  to  separate  the  amplitudes  of  the  fundamental  disturbance  frequency  and  its  higher  har¬ 
monics,  as  mentioned  earlier.  Figures  4.22  through  4.25  show  the  mean  flow  and  the  amplitude 
distribution  over  the  radius  for  q  of  the  fundamental,  the  first,  and  the  second  higher  harmonic 
at  the  cross  section  a;  =  3.5.  As  these  figures  show,  the  results  compare  very  well,  in  spite  of 
the  different  disturbance  methods  used  in  the  experiment  and  the  simulation.  At  this  point  it  is 
unclear  why  the  maximum  at  the  axis  of  symmetry  for  the  first  harmonic  shows  appears  to  be 
lower  for  the  numerical  simulation.  A  possible  explanation  could  be  the  low  resolution  near  the 
axis  due  to  the  grid  stretching  or  insufficient  data  points  in  time  for  the  Fourier  transformation. 

In  summary,  it  was  found  that  with  the  compressible  code  it  is  possible  to  accurately  calcu¬ 
late  axisymmetric  and  three-dimensional  flow  fields  for  a  subsonic  Mach  number  of  Mqo  =  0.2. 
Exceptionally  good  agreement  was  found  with  the  water  channel  experiments  by  Siegel  (1996) 
and  the  incompressible  simulations  by  [Schwarz  (1996)].  Further,  it  was  found  that  the  com¬ 
pressible  simulation  can  be  a  valuable  tool  for  the  calculation  of  incompressible  wahe  flow  fields. 
As  mentioned  earlier,  all  simulations  were  performed  on  a  single  processor  of  a  Silicon  Graphics 
Inc.  Power  Challenge  L  with  4  R8000  CPUs  (90  MHz,  4MB  cache  each)  and  512  MBytes  of 


mam  memory. 
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Figure  4.21:  Isolines  of  the  fundamental  of  first  azimuthal  mode  of  azimuthal  velocity  for  three 
difllerent  spatial  resolutions  for  Mqo  =  0.2  and  Beo  =  1,000  [cases  T3  (top)  through  T5  (bot- 


o  experiment 
-  simulation 


Figure  4.22:  Mean  fiow  profiles  for  Mqo  =  0.2  and  ife/j  =  1,000  at  a;  =  3.5. 


Figure  4.23:  Normalized  amplitude  distribution  of  fundamental  disturbance  frequency  for  M^o  = 
0.2  and  Bed  =  1,000  at  x  =  3.5. 


y/D 


Figure  4.24:  Normalized  amplitude  distribution  of  first  higher  harmonic  of  disturbance  frequency 
for  Moo  =  0.2  and  ifep  =  1, 000  at  x  =  3.5. 


Figure  4.25:  Normalized  amplitude  distribution  of  second  higher  harmonic  of  distmrbance  fre 
quency  for  M^o  =  0.2  and  Beo  =  1,000  at  a;  =  3.5. 
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4.3  Large-Eddy  Simulation 

A  preliminary  validation  for  the  LES  code  was  performed  for  a  free  stream  Mach  number  of 
Moo  =  0.2  and  a  global  Reynolds  number  of  Beo  =  2,000.  As  before,  the  results  were  directly 
compared  qualitatively  to  experimental  measurements  obtained  in  the  water  channel  of  the 
Hydrodynamics  Laboratory  of  the  Aerospace  and  Mechanical  Engineering  Department  at  the 
University  of  Arizona  by  [Siegel  (1996)].  The  experiments  have  shown  that  the  flow  is  fully 
turbulent  and  exhibits  large  coherent  structures,  which  are  very  similar  to  the  laminar  structures 
that  appeared  at  ife/j  =  1,000. 

For  the  numerical  simulation,  at  flrst  the  flow  fleld  was  calculated  with  DNS  in  exactly  the 
same  way  as  the  flow  at  ife^)  =  1, 000.  Dining  the  initial  stage  a  periodic  disturbance  in  the  flrst 
azimuthal  mode  was  generated  through  the  blowing  and  suction  slot.  After  flve  time  periods 
the  disturbance  was  ramped  down  to  zero.  Figure  4.26  shows  the  time  response  of  the  flow  at 
X  =  0.25  and  r  =  0.25  in  form  of  the  total  azimuthal  vorticity  for  the  flrst  flve  periods  in  time. 
As  for  the  lower  Reynolds  number  case,  the  disturbance  grows  exponentially  in  time. 


Figure  4.26:  Initital  time  response  of  total  azimuthal  vorticity  to  a  single  pulse  disturbance  of 
Fourier  mode  fc  =  1  at  r=0.25  and  z=  0.25. 

In  contrast  to  the  lower  Reynolds  number  case,  however,  the  structures  that  appear  in  the 
free  shear  layer  became  too  small  for  the  spatial  resolution  of  the  DNS.  This  resulted  in  high 
local  gradients.  Shortly  after  t  =  72  the  calculation  was  underresolved  and  terminated.  At  this 
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point  the  simulation  was  continued  with  the  LES  code,  using  the  same  spatial  resolution  as 
before  for  the  DNS.  Figure  4.27  shows  the  continued  response  of  the  flow  at  the  same  location 
as  before.  As  for  the  laminar  case,  the  response  remains  highly  unsteady  and  non-periodic  even 
though  the  flow  is  not  disturbed  any  further. 


Figure  4.27:  Time  response  of  total  azimuthal  vorticity  to  a  single  pulse  disturbance  of  Fourier 
mode  A:  =  1  at  r=0.25  and  z=  0.25. 

At  this  point,  large  coherent  structures  appear  as  shown,  for  example,  in  flgme  4.28  in  the 
form  of  isolines  of  instantaneous  total  vorticity.  The  structures  look  very  similar  to  the  laminar 
structures  for  the  lower  Reynolds  number  case.  In  addition,  flow  visualization  via  particles  shows 
that  the  larger  structmes  are  similar  to  the  structures  that  were  observed  in  the  water  channel 
experiments  (see  figmre  4.29). 

More  code  validation  for  the  LES  code  is  needed  at  this  point,  especially  for  supersonic  flows. 
So  far,  flow  flelds  at  a  supersonic  Mach  number  of  Moo  =  2.46  and  global  Reynolds  numbers  up 
to  Ben  =  400, 000  have  been  calculated  (see  chapter  7).  In  order  to  perform  a  direct  comparison 
with  experiments,  however,  flows  at  higher  Reynolds  numbers  {Esd  =  0(1,000,000))  need  to 
be  calculated. 

Therefore,  in  the  present  work  only  preliminary  results  are  presented  for  the  LES  calculations 
of  the  higher  Reynolds  number  (turbulent)  supersonic  flow  flelds. 
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Figure  4.29:  Flow  visualization  using  particles  for  Moo  —  0-2  and  ifep  =  2, 000;  simulation  (top) 
and  water  channel  experiment  by  Siegel  (1994)  (bottom). 


Chapter  5 


Results 


Three-dimensional  unsteady  flow  fields  of  the  wake  of  an  axsymmetric  blufiF  body  with  a  blunt 
base  have  been  calculated  for  several  combinations  of  free  stream  Mach  number  and  global 
Reynolds  number.  The  main  goal  for  the  present  calculations  was  to  investigate  the  existence 
of  large  structures  and  their  influence  and  significance  on  the  characteristics  of  the  global  flow 
field.  As  a  first  step,  Direct  Numerical  Simulations  (DNS)  were  performed  to  investigate  the  time 
dependent  behavior  of  laminar  flows  in  order  to  determine  the  existence  of  absolute  instabilities 
and  investigate  the  evolution  of  large  structures  in  laminar  wake  flows  and  their  effect  on  the 
global  flow  field.  Secondly,  although  preliminary,  the  calculations  were  extended  to  turbulent 
flows  by  employing  Large-Eddy  Simulations  (LES)  in  order  to  investigate  if  absolute  instabilities 
also  exist  and  if  large  coherent  structures  are  also  present  in  turbulent  supersonic  wake  flows. 

5.1  Direct  Numerical  Simulations 

As  pointed  out  earlier,  the  incompressible  wake  of  axisymmetric  bodies  at  a  global  Reynolds 
number  of  ifep  =  1, 000  is  subject  to  an  absolute  instability  in  the  first  helical  Fourier  mode.  As  a 
result,  the  flow  is  highly  unsteady  and  dominated  by  large  structures.  However,  the  flow  remains 
laminar  at  all  times  and  does  not  undergo  transition  to  turbulence.  Prom  previous  investigations 
[see,  for  example,  [Schwarz  (1996)]  and  Siegel  (1994)]  there  is  a  basis  of  knowledge  about  the 
evolution  of  large  structures  in  the  incompressible  wake.  For  that  reason,  the  low  Mach  number 
subsonic  flow  has  been  used  for  code  validation  and  is  also  used  here  for  comparison  of  the 
evolution  of  the  structures  in  supersonic  flows.  In  the  following  sections  results  from  DNS  are 
presented  for  three  different  free  stream  Mach  numbers,  M^o  =  0.2,  Mqq  =  1.2,  and  M^o  =  2.46. 

5.1.1  Subsonic  Flow 
a)  Steady  Calculations 

The  subsonic  flow  field,  that  has  been  evaluated  for  validation  of  the  DNS  code  (cases  A1 
and  T3  through  T5)  is  presented  here  to  provide  a  basis  for  comparison  of  the  results  with 
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supersonic  wake  flows.  First,  additional  information  on  the  steady  axisymmetric  flow  field  is 
provided  in  form  of  isolines  of  axial  velocity,  radial  velocity,  pressure,  density,  and  temperature. 
These  results  are  shown  in  figures  5.1  through  5.5. 


Figure  5.1:  Isolines  of  axial  velocity  for  Moo  =  0.2  and  Ren  =  1,000  (axisymmetric  unforced 
calculation). 


Figure  5.2:  Isolines  of  radial  velocity  for  Mo©  =  0.2  and  Iho  =  1,000  (axisymmetric  unforced 
calculation). 


Figure  5.3:  Isolines  of  pressure  for  Moo  =  0.2  and  ifex?  =  1,000  (axisymmetric  unforced  calcula¬ 
tion). 

The  most  important  feature  at  this  point  can  be  seen  in  the  pressure  distribution  (figure  5.3. 
The  graph  clearly  shows  the  strong  pressure  drop  immediately  downstream  of  the  base,  which 
is  responsible  for  the  base  drag.  However,  as  will  be  seen  in  the  following  sections,  the  pressure 
drop  is  much  more  pronounced  in  supersonic  flow  fields.  In  addition,  the  plot  reveals  an  elevated 
pressure  at  the  base  near  the  axis  of  symmetry.  This  can  be  explained  by  the  stagnation  point 
flow  behavior  of  the  recirculating  region.  This  results  in  the  typical  base  pressure  distribution 
as  shown  in  figure  5.19. 


b)  Unsteady  Calculations 
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Figure  5.4:  Isolines  of  density  for  Moo  =  0.2  and  Ren  =  1, 000  (axisymmetric  unforced  calcula¬ 
tion). 


Figiure  5.5:  Isolines  of  temperature  for  Moo  =  0.2  and  Bbd  =  1,000  (axisymmetric  unforced 
calculation). 

As  mentioned  in  the  discussion  concerning  code  validation  in  section  4.3,  for  the  DNS  of  the 
three-dimensional  unsteady  wake  the  flow  was  disturbed  by  a  single  pulse  disturbance  locally 
within  the  recirculating  region  at  the  location 

®dist  =  0-25  and  r^ist  =  0.25,  (5.1) 

as  described  in  section  3.7.  The  time  dependent  response  of  the  flow  field  at  the  location 
of  the  disturbance  is  shown  in  figme  4.13.  After  a  period  of  time,  the  signal  deviates  from 
the  axisymmetric  state  and  starts  oscillating  in  a  non-periodic  fashion.  At  this  point,  large 
structures  are  present  in  the  flow. 

Four  difierent  methods  are  used  here  for  the  identification  of  vortical  structures  in  the  flow 
field: 

•  Isolines  of  instantaneous  total  vorticity  shown  in  difierent  planes,  the  x-y  plane  and  the  x-z 
plane  through  the  axis  of  symmetry,  and  the  y-z  plane  at  4  difierent  x  locations  (x  =  1.0, 
X  =  2.0,  X  =  3.0,  and  x  =  4.0). 

•  Visualization  of  the  instantaneous  flow  field  using  particles. 

•  Isosurfaces  of  instantaneous  total  vorticity,  axial  vorticity,  and  azimuthal  vorticity. 

•  Isosurfaces  of  instantaneous  pressure  deviation  from  the  time-averaged  flow  field. 

Figure  5.6  shows  isolines  of  instantaneous  total  vorticity  in  two  planes  perpendicular  to  each 
other.  The  figure  reveals  that  the  flow  remains  symmetric  in  the  second  plane.  This  plane 
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Figure  5.6:  Instantaneous  total  vorticity  for  Ren  =  1,000  and  M^o  =  0.2.  x-z  plane  (top)  and 
x-y  plane  (bottom). 

of  symmetry  has  also  been  observed  in  the  incompressible  calculations  and  in  water  channel 
experiments.  Therefore,  in  order  to  save  computer  resources  a  plane  symmetry  was  prescribed 
in  all  following  calculations. 

Strong  vortical  structures  appear  which  originate  near  the  corner  of  the  base  and  are  con- 
vected  downstream.  These  structures  have  also  been  observed  in  water  channel  experiments. 
Figure  5.7  shows  isolines  of  the  total  vorticity  at  axial  cross-sections  at  four  different  x  locations. 
The  figures  clearly  indicate  two  strong  maxima  of  total  vorticity  at  every  cross-section.  However, 
some  cross-sections  exhibit  more  than  two  local  maxima  (for  example  at  a:  =  3).  The  structure 
of  the  axial  vortices  will  become  clearer  when  looking  only  at  the  axial  vorticity,  which  is  shown 
in  the  next  series  of  figures. 

Figures  5.8  through  5.10  show  isosmfaces  of  total  vorticity,  azimuthal  vorticity  and  total 
axial  vorticity,  respectively.  Comparing  the  total  vorticity  with  the  azimuthal  vorticity,  which 
are  both  shown  at  the  same  level  of  lw|  =  3.0,  reveals  that  the  ring-shaped  structure  consists 
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Figure  5.7:  Instantaneaus  total  vorticity  for  ifep  =  1,000  and  Moo  =  0.2,  y-z  planes  at  a;  =  1 
(top  left),  X  =  2  (top  right),  a:  =  3  (bottom  left),  and  a;  =  4  (bottom  right). 

mainly  of  azimuthal  vorticity.  This  structure  is  followed  by  two  structures  of  axial  vorticity 
(shown  here  at  \ux\  =  1.0),  which  are  counter-rotating.  The  connection  between  the  two  types 
of  structures  becomes  clearer  when  looking  at  figure  5.11.  Here  the  isosurface  of  instantaneous 
pressure  deviation  from  the  time-averaged  flow  field  is  shown  (for  Ap  =  —0.02).  The  figure 
shows,  that  two  ring-shaped  structures  are  connected  by  two  axial  structures. 


CHAPTERS.  R 


75 


Figure  5.9:  Isosurfaces  of  instantaneous  axial  vorticity  (Iwxl  =  l-O)  for  .^oo  =  0.2,  Rbd  =  1,000. 

For  a  direct  comparison  with  flow  visualization  in  water  channel  experiments  that  were 
carried  out  by  Siegel  (1994),  particles  have  been  introduced  into  the  flow  during  the  numerical 
simulations.  The  locations  of  the  particles  are  updated  at  every  time  step  during  the  numerical 
simulation.  The  release  of  new  particles  is  triggered  at  a  specified  number  of  time  steps  (typically 
every  16  to  64  steps).  The  particles  are  released  at  a  ring  near  the  corner  of  the  base,  such  that 
they  enter  the  recirculation  region  close  to  the  center  of  the  shear  layer.  This  resembles  the 
introduction  of  dye  at  the  corner  of  the  base  in  the  water  channel  experiments. 

In  figure  4.14  the  flow  field  is  shown  for  two  side  views  perpendicular  to  each  other.  As  before, 
the  figures  show  the  symmetry  to  the  x-y  plane  and  the  appearance  of  a  ring-like  structure 
followed  by  two  axial  vortices.  In  figure  4.15  particles  are  shown  at  x  locations  x  =  1.0,  x  =  2.0, 
X  =  3.0,  and  x  =  4.0.  This  method  is  similar  to  a  laser  sheet  illumination  of  fluorescent  dye 
in  a  water  channel  experiment.  Again  the  figmres  reveal  the  symmetry  to  the  x-y  plane  and 
clearly  shows  the  two  rollups  which  result  from  a  superposition  of  a  A:  =  1  and  a  A:  =  —  1 
helical  disturbance  mode.  Similar  structures  have  been  observed  in  water  channel  experiments 
when  using  a  laser  sheet  illumination  technique  at  cross  sections  in  the  y-z  plane  (compare,  for 
example,  figure  5.17). 

The  three-dimensional  time  dependent  results  show  that  there  is  a  strong  unsteady  part  in 
the  flow.  Large  structures  are  present.  In  order  to  investigate  the  influence  of  the  structures  on 
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Figure  5.10;  Isosurfaces  of  instantaneous  aaimuthal  vorticity  =  3.0)  for  M^o  =  0.2,  Iho  = 
1,000. 


Figure  5.11:  Isosurfaces  of  instantaneous  pressure  deviation  from  the  time-averaged  mean  flow 
{Ap  =  -0.02)  for  Moo  =  0.2,  Iko  =  1,000. 


the  global  flow  fleld,  a  time  average  is  performed  over  a  time  period  from  t  =  80  to  f  =  150. 
Figures  5.12  through  5.18  show  isolines  of  axial  velocity,  radial  velocity,  pressmre,  density,  tem¬ 
perature,  azimuthal  vorticity  and  stream  lines  for  the  time-averaged  flow  fleld,  respectively.  The 
strongest  changes  compared  to  the  axisymmetric  steady  flow  are  noticeable  in  the  distribution 
of  the  pressure  and  the  density.  The  streamlines  reveal  that  the  recirculation  region  for  the 
time-averaged  flow  field  is  only  1.8  diameters  long,  compared  to  the  recirculation  length  of  2.1 
diameters  for  the  steady  flow. 

The  pressure  distribution  along  the  base  has  become  almost  constant,  in  contrast  to  the 
steady  calculation  which  exhibits  a  strong  maximum  at  the  axis  of  symmetry  (see  figure  5.19). 
This  phenomenon  is  even  more  pronounced  in  the  case  of  supersonic  flows,  as  will  be  shown  in 
the  following  sections. 
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Figure  5.12:  Isolines  of  axial  velocity  for  the  time-averaged  flow  of  M^o  =  0.2  and  Ren  =  1,000. 


Figure  5.13:  Isolines  of  radial  velocity  for  the  time-averaged  flow  of  Mqo  =  0.2  and  ifejj  =  1, 000 


Figure  5.14:  Isolines  of  pressure  for  the  time-averaged  flow  of  Moo  =  0.2  and  Reu  =  1,000 


Figure  5.16:  Isolines  of  temperature  for  the  time-averaged  flow  of  Mqo  =  0.2  and  Rej)  =  1,000 
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5.1.2  Supersonic  Flow  at  Mach  1.2 

For  the  first  simulation  of  a  supersonic  wake  the  free  stream  Mach  number  was  Moo  =  1-2  and 
the  global  Reynolds  number  was  Ren  =  4,000  (cases  A6  and  T6). 

a)  Steady  Calculations 

As  for  the  subsonic  case,  at  first  a  steady  axisymmetric  fiow  was  calculated.  The  resulting 
flow  field  is  shown  in  figures  5.20  through  5.26  in  form  of  isolines  of  axial  velocity,  radial  velocity, 
streamlines,  pressure,  density,  temperatme,  and  azimuthal  vorticity,  respectively.  The  figures 
indicate  that  the  features  of  the  supersonic  flow  are  different  from  its  subsonic  counterpart. 


Figure  5.20:  Isolines  of  axial  velocity  for  Moo  —  1-2  and  Rbd  =  4,000  (axisymmetric  steady 
calculation). 


Figure  5.21:  Isolines  of  radial  velocity  for  Moo  =  1-2  and  Red  =  4,000  (axisymmetric  steady 
calculation). 

The  figures  for  the  axial  velocity  do  not  show  significant  differences  from  the  subsonic  flow. 
The  approaching  boundary  layer  was  thinner  at  the  corner  for  the  supersonic  flow,  as  can  be 
seen  in  figure  5.20.  This  is  determined  by  the  prescribed  inflow  condition.  The  radial  velocity, 
however,  looks  different  for  the  different  Mach  numbers.  Here,  the  area  of  negative  velocity 
above  the  free  shear  layer  extends  farther  upstream  than  for  the  subsonic  calculation.  It  also 
reveals  the  existence  of  expansion  waves.  The  streamlines,  shown  in  figure  5.22,  have  very  similar 
shapes  to  the  ones  of  the  subsonic  flow,  with  the  exception  that  the  curvature  of  the  separating 
streamline  is  much  weaker  for  the  supersonic  case. 

The  largest  difference  can  be  observed  in  the  pressure  distribution.  Figure  5.23  clearly  shows 
the  regions  of  expansion  and  of  weak  re-compression  which,  of  course,  are  not  present  in  the 
subsonic  flow.  In  addition,  the  change  in  pressme  is  much  stronger  for  the  supersonic  flow.  This 
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Figure  5.22:  Steady  streamlines  for  Moo  =  1-2  and  Rej)  =  4,000  (axisymmetric  steady  calcula¬ 
tion). 


Figure  5.23:  Isolines  of  pressure  for  Mqo  =  1.2  and  Reu  =  4, 000  (axisymmetric  steady  calcula¬ 
tion). 


results  in  a  much  higher  base  pressure  drag,  which  agrees  with  observations  from  experiments 
(see,  for  example,  [Rollstin  (1987)]).  Similar  features  as  in  the  pressure  can  be  seen  in  the  density 
and  temperatme  distribution,  which  both  show  a  much  stronger  variation  for  the  supersonic  flow 
than  for  the  subsonic  flow. 

The  isolines  of  azimuthal  vorticity  look  very  similar  to  those  of  the  subsonic  flow.  However, 
the  absolute  values  of  the  vorticity  are  much  higher  for  the  present  case. 

b)  Unsteady  Calculations 

Similar  to  the  subsonic  calculation,  as  discussed  in  the  preceding  section,  the  flow  was  dis¬ 
turbed  locally  with  a  single  pulse  disturbance.  The  location  of  the  disturbance  introduction  was 
at  r  =  0.25  and  x  =  0.25.  Figme  5.27  shows  the  time  signal  of  the  radial  momentum  at  the 
location  of  the  disturbance  introduction  for  Fourier  modes  k  =  0,  k  =  1,  and  k  =  2.  The  signal 
exhibits  a  behavior  similar  to  that  observed  for  the  subsonic  calculation.  After  some  time  the 
signal  deviates  from  the  axisymmetric  state  and  starts  oscillating  in  a  non-periodic  fashion. 

On  a  logarithmic  scale,  which  is  shown  in  figure  5.28,  the  signal  of  the  radial  momentum  at 
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Figure  5.24:  Isolines  of  density  for  Mqo  =  1.2  and  Rbd  =  4,000  (axisymmetric  steady  calcula¬ 
tion). 


Figure  5.25:  Isolines  of  temperature  for  Mqo  =  1.2  and  Rep  =  4,000  (axisymmetric  steady 
calculation). 


r  =  0.25  and  x  =  0.25  exhibits  an  exponential  deviation  from  the  axisymmetric  state  until  a  state 
of  nonlinear  saturation  is  reached.  This  is  a  clear  indication  that  an  absolute  instability  exists. 
In  the  following,  the  response  of  the  flow  is  further  studied  in  the  same  fashion  as  discussed  for 
the  subsonic  flow  fleld.  Figmres  5.29  through  5.31  show  isolines  of  instantaneous  total  vorticity 
for  difierent  time  instances.  This  series  reveals  the  development  of  vortical  structmes  which 
originate  near  the  corner  of  the  base  and  increase  in  strength  while  being  convected  downstream. 
The  structmes  appear  to  be  similar  to  the  ones  observed  in  the  subsonic  case,  except  that  they 
have  a  much  higher  concentration  of  vorticity.  In  addition,  a  layer  of  high  vorticity  develops 
near  the  center  line  farther  upstream  (see  flgure  5.31). 

As  for  the  subsonic  case,  the  figmes  reveal  the  existence  of  a  plane  of  symmetry  as  shown  in 
flgmres  5.32  through  5.34  in  which  isolines  of  instantaneous  total  vorticity  in  the  x-y  plane  are 


Figure  5.26:  Isolines  of  azimuthal  vorticity  field  for  Mqo  =  1.2  and  Ben  =  4,000  (axisymmetric 
steady  calculation). 
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Figure  5.27:  Time  response  of  radial  momentum  to  a  single  pulse  of  Fourier  mode  fc  =  1  at 
r  =  0.25  and  x  =  0.25  at  Moo  =  1.2  and  ifejp  =  4, 000. 


plotted  at  four  different  x  locations.  These  plots  reveal  that  in  addition  to  the  high  concentration 
of  vorticity  at  two  locations  in  a  cross  section,  locations  of  high  concentration  of  vorticity  are  also 
present  closer  to  the  axis.  These  phenomena  were  also  observed  for  the  subsonic  case.  The  three- 
dimensional  shape  of  the  vortical  structures  can  be  seen  more  clearly  in  figures  5.35  through  5.38, 
where  isosurfaces  of  instantaneous  total  vorticity,  azimuthal  vorticity,  axial  vorticity  and  pressure 
deviation  from  the  time-averaged  mean  flow  are  plotted,  respectively.  As  for  the  subsonic  case, 
the  total  vorticity  and  the  azimuthal  vorticity  exhibit  ring-shaped  structures,  which  are  followed 
by  a  pair  of  axial  vortices  of  lesser  strength  (see  figure  5.36). 

Flow  visualization  by  particle  introduction  as  shown  in  figures  5.39  through  5.40,  also  reveals 
the  structures.  As  before,  the  particles  were  introduced  at  a  ring  near  the  corner  of  the  base. 
However,  the  structures  are  not  as  clearly  noticeable  as  in  the  subsonic  case.  To  date  there 
have  been  no  visualization  experiments  for  flows  at  this  Mach  number.  Therefore,  no  direct 
comparison  is  possible  to  experimental  observation. 
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Figure  5.28:  Time  response  of  radial  momentum  to  a  single  pulse  of  Fourier  mode  /c  =  1  at 
r  =  0.25  and  x  =  0.25  on  a  logarithmic  scale  at  Mqo  =  1.2  and  Ren  =  4, 000. 


Figure  5.29:  Isolines  of  instantaneous  total  vo^ticity  for  Ben  =  4,000  and  Moo  =  1-2,  at  0  =  0 
(top)  and  6  =  7r/2  (bottom)  at  t  =  121.34. 
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Figure  5.36:  Isosurfaces  of  instantaneous  axial  vorticity  (|n)a;|  =  2)  for  M^o  =  1.2,  =  4,000 

at  t  =  135.09. 


Figure  5.37;  Isosurfaces  of  instantaneous  azimuthal  vorticity  [\uJe\  =  6)  for  =  1.2,  ife/j  = 
4,000  at  t  =  135.09. 


Figure  5.38:  Isosurfaces  of  instantaneous  pressure  deviation  from  the  time-averaged  mean  flow 
(Ap  =  —0.02)  for  Moo  =  1.2,ife£)  =  4,000  at  t  =  135.09. 
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Figure  5.39:  Flow  visualization  of  instantaneous  flow  field  via  particles  for  Moo  =  1'2  and 
Ben  =  4,000,  viewed  from  6  =  7r/2  (top)  and  from  0  =  0  (bottom)  (particles  introduced  at 
a:  =  0.025,r  =  0.475). 


Figure  5.40:  Flow  visualization  of  instantaneous  flow  field  via  particles  for  Mqo  =  1.2  and 
Red  =  4, 000,  sheets  of  thickness  0.01  at  0  =  0  (top)  and  6  =  7r/2  (bottom),  (particles  introduced 
at  x  =  0.025,r  =  0.475). 


The  subsonic  calculations  have  shown  that  the  unsteady  structures  can  have  a  strong  in¬ 
fluence  on  the  global  time-averaged  flow.  This  influence  was  strongest  in  the  time-averaged 
pressure  and  the  density.  For  the  case  of  M^o  =  1.2,  isolines  of  the  time-averaged  pressure  and 
density  are  shown  in  figures  5.42  and  5.43,  respectively.  They  indicate  that  the  distribution  of 
the  pressure  and  the  density  along  the  base  is  almost  constant.  This  finding  is  in  good  agreement 
with  experimental  observations  of  turbulent  supersonic  axisymmetric  wake  flows.  The  difference 
between  the  steady  axisymmetric  and  the  time-averaged  base  pressure  is  again  obvious  from  fig¬ 
ure  5.44.  Here  the  pressure  is  normalized  by  the  free  stream  static  pressure.  The  plot  shows 
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Figure  5.41:  Flow  visualization  of  instantaneous  flow  fleld  via  particles  for  Mqo  =  1.2  and 
Rejj  =  4,000,  cross-sectional  cuts  of  thickness  0.01  at  x  =  1.0,  x  =  2.0,  x  =  3.0,  and  x  =  4.0 
(particles  introduced  at  x  =  0.025,  r  =  0.475). 


Figure  5.42:  Isolines  of  pressiure  for  a  time-averaged  flow  at  Moo  =  1-2  and  ife/j  =  4,000. 

that,  in  addition  to  a  change  in  shape  of  the  pressure  distribution,  the  average  pressure  is  lower 
than  the  average  of  the  steady  axisymmetric  pressiure.  This  results  in  a  larger  base  drag,  which 
is  caused  solely  by  the  unsteady  structures. 

In  summary,  the  supersonic  flow  for  Moo  =  1.2  and  Rbd  =  4,000  shows  a  behavior  that 
is  somewhat  similar  to  that  presented  previously  for  the  subsonic  case.  The  flow  exhibits  an 
absolute  instability,  which  results  in  the  formation  of  unsteady  structures.  These  structures  have 
a  strong  influence  on  the  pressure  distribution  of  the  time-averaged  flow,  resulting  in  a  constant 
base  pressure  distribution  as  found  in  experiments,  and  a  higher  base  drag  than  for  the  steady 
axisymmetric  calculation.  The  absolute  instability  causes  an  exponentially  growing  deviation 
from  the  axisymmetric  flow.  After  a  period  of  time  the  deviation  reaches  a  nonlinear  saturation 
and  the  flow  field  fluctuates  in  a  non  periodic  way.  The  unsteady  flow  remains  laminar  and  does 
not  show  any  sign  of  a  transition  to  turbulence. 
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Figure  5.43:  Isolines  of  density  for  a  time-averaged  flow  at  Mq©  =  1-2  and  Her)  =  4,000. 


r 


Figure  5.44:  Base  pressure  distributions  for  M©©  =  1.2  and  Red  =  4,000. 

5.1.3  Supersonic  Flow  at  Mach  2.46 

The  flow  field  of  case  A4,  which  was  calculated  for  the  validation  of  the  DNS  code  for  axisym- 
metric  flow  calculations  (see  section  4.1.2),  was  used  as  initial  condition  for  further  investigations 
of  the  time  dependent  behavior  of  the  three-dimensional  flow. 

a)  Steady  Calculations 

In  addition  to  the  results  presented  for  the  steady  axisymmetric  flow  field  in  chapter  4, 
figures  5.45  and  5.46  show  the  streamlines  and  isolines  of  vorticity,  respectively.  As  pointed  out 
in  the  previous  section  for  the  flow  with  free  stream  Mach  number  M©©  =  1.2,  the  plot  of  the 
streamlines  reveals  that  the  curvature  of  the  streamline  dividing  the  recirculating  region  from 
the  outer  region  is  much  smaller  for  supersonic  flow.  For  the  present  case  at  M©©  =  2.46,  the 
dividing  streamline  is  almost  a  straight  line. 

Also,  as  already  mentioned  in  the  previous  section,  the  pressure  distributipn  is  significantly 
difierent  for  the  subsonic  and  supersonic  case.  This  difference  becomes  even  more  pronounced 
for  the  higher  Mach  number  (comparison  of  figure  4.9  with  figures  5.23  and  5.3).  The  overall 


CHAPTERS.  RESULTS 


94 


variation  of  the  pressure  is  much  larger  for  this  case  than  for  the  previously  shown  results.  The 
same  effect  can  be  seen  in  the  base  pressure  variation,  which  results  in  a  higher  base  drag  for 
higher  free  stream  Mach  numbers.  This  is  indicated  by  the  base  pressure  distributions  shown  in 
figures  5.74,  5.44,  and  5.19.  The  azimuthal  vorticity,  however,  shows  very  similar  characteristics 


Figure  5.45:  Steady  streamlines  for  Moo  =  2.46  and  ifep  =  30,000  (axisymmetric  steady  calcu¬ 
lation)  . 
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b)  Unsteady  Calculations 

In  calculations  for  a  free  stream  Mach  number  of  M^o  =  2.46  [Tourbier  &  Fasel  (1994)],  no 
absolute  instability  could  be  observed  up  to  a  global  Reynolds  number  of  Red  =  3,000.  Even 
when  increasing  the  Reynolds  number  up  to  ifejo  =  25,000  in  subsequent  calculations,  there 
was  no  indication  of  an  absolute  instability.  For  all  cases  with  Reynolds  numbers  lower  than 
Red  =  30, 000  a  generated  pulse  distmbance  always  decayed  in  time  and  the  flow  fleld  eventually 
relaxed  to  the  original  axisymmetric  steady  state.  The  unsteady  flow  behavior  in  all  these  cases 
was  similar  to  that  presented  by  Tourbier  &  Fasel  (1994).  The  main  difierence  was  that  the 
temporal  decay  rates  of  the  disturbances  decreased  for  increasing  Reynolds  numbers.  This  is 
an  indication  that  the  flow  becomes  less  stable  when  the  Reynolds  number  is  increased.  Here, 
results  of  a  DNS  are  presented  a  global  Reynolds  number  Red  =  30,000. 


Figure  5.47:  Time  response  of  radial  momentum  to  a  single  pulse  of  Fourier  mode  fc  =  1  at 
r  =  0.25  and  x  =  0.25  for  Moo  =  2.46  and  Red  =  30,000  (unfiltered  DNS). 

The  axisymmetric  flow  was  disturbed  with  a  single  pulse  disturbance  at  a:  =  0.25  and  r  = 
0.25.  The  time  signal  at  the  location  of  the  disturbance  introduction  is  shown  in  figure  5.47.  As 
for  both  the  subsonic  case  and  the  supersonic  case  with  Moo  =  1-2,  the  azimuthal  decomposition 
of  the  time  signal  exhibits  an  exponential  deviation  from  the  axisymmetric  state.  Figure  5.47 
shows  that  the  Fomier  mode  fc  =  1  is  not  the  most  amplified  for  this  flow;  k  =  2,  k  =  3, 
and  k  =  4  show  stronger  amplification  rates.  When  the  first  four  Fourier  modes  pass  a  certain 
threshold,  the  higher  Fourier  modes  (A:  =  5,  fc  =  6,  A:  =  7,  and  k  =  8)  also  exhibit  strong 
amplification. 

In  contrast  to  the  two  cases  discussed  previously,  local  velocity  gradients  near  the  stagnation 
point  (x  w  2.8)  kept  increasing  over  time  to  very  high  levels.  Eventually,  the  simulation  became 
underresolved  and  terminated  at  about  t  =  125.  Figures  5.48  and  5.49  show  instantaneous 
isolines  of  total  vorticity  shortly  before  the  simulation  terminated.  At  this  point  the  flow  field 
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Figure  5,48:  Isolines  of  instantaneous  total  vorticity  for  =  30,000  and  Mo©  ==  2.46.  x-z 
plane  (top)  and  x-y  plane  (bottom)  (unfiltered  DNS). 

exhibits  large  and  small  scale  structures  (see  figure  5.48).  The  simulation  is  no  longer  able  to 
handle  the  large  gradients  in  the  flow  field,  which  causes  numerical  instabilities  and  finally  a 
termination  of  the  calculation.  Also,  as  can  be  seen  in  figure  5,47,  all  azimuthal  Fourier  modes 
fill  up  rapidly  and  even  the  higher  modes  grow  up  to  the  same  order  of  magnitude.  All  of  this  is 
an  indication  that  the  spatial  resolution  of  the  simulation  is  inadequate  for  the  occurring  small 
scale  structures. 

Because  of  the  appearance  of  small  scale  structures  that  could  not  be  resolved  with  a  feasible 
number  of  grid  points  and  Fourier  modes  (for  a  reasonable  computation  time  on  the  computer 
that  was  used)  the  flow  field  was  filtered  in  the  axial  amd  radial  directions  at  every  time  step  of  the 
Runge-Kutta  integration  scheme.  For  this  procedure  a  sixth-order  accurate  compact  difference 
filter  was  applied,  which  is  described  in  Lele  (1992).  This  filter  has  been  used  successfully 
for  Large-Eddy  Simulations  of  incompressible  boundary  layer  flows  [Bachman  (1996)]  and  has 
also  been  applied  for  LES  in  the  present  work  (see  section  3.6).  In  addition,  the  flow  field 
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Figure  5.49:  Isolines  of  instantaneous  total  vorticity  for  ihjj  =230,000  and  M^o  =  2.46,  y-z 
planes  at  a:  =  1  (top  left),  x  =  2  (top  right),  x  =  3  (bottom  left),  and  x  =  4  (bottom  right) 
(unfiltered  DNS). 


was  also  filtered  at  every  hundredth  time  step  with  an  explicit  second  order  central  difference 
filter,  because  the  dissipation  provided  by  the  sixth-order  compact  filter  was  not  sufficient  to 
completely  suppress  the  growth  of  small  scales. 

The  time  signal  at  the  location  of  the  disturbance  generation  after  the  introduction  of  the 
spatial  filtering  is  shown  in  figme  5.50.  Even  with  the  added  dissipation  of  the  filtering  the  flow 
field  still  exhibits  non-periodic  fluctuations  in  all  Fourier  modes  that  do  not  decay.  In  figures 
5.51,  5.52,  and  5.53  isolines  of  total  vorticity  are  plotted  in  the  x-y  plane  and  the  x-z  plane 
for  the  spatially  filtered  calculations  for  three  different  instances  in  time.  As  for  the  previous 
calculations,  the  figures  vortical  structures  appear  in  the  free  shear  layer  which  originate  near 
the  corner  of  the  base  and  are  convected  downstream.  In  contrast  to  the  previous  cases,  the 
structures  are  more  confined  to  the  shear  layer  region.  This  is  confirmed  by  figures  5.54,  5.55, 
and  5.56  which  show  isolines  of  total  vorticity  in  the  y-z  plane  at  four  different  x  locations. 
The  structures  that  are  similar  to  the  ones  observed  for  the  lower  Mach  number  cases  are  now 
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Figure  5.50:  Time  signal  of  radial  momentum  at  r  =  0.25  and  x  =  0.25  for  Moo  =  2.46  and 
/fee  =  30,000  (filtered  DNS). 

located  entirely  within  the  top  shear  layer.  The  three-dimensional  structures  can  be  seen  more 
clearly  in  the  three-dimensional  isosurfaces  in  figures  5.57  through  5.68.  The  isosurface  for  the 
instantaneous  pressure  deviation  from  the  time-averaged  flow  is  shown  at  Ap  =  —0.02  which 
represents  the  same  level  as  for  the  previous  cases.  It  appears  that  there  are  still  two  axial 
vortices  present.  The  vorticity  structures  indicate  that  for  this  case  the  level  of  vorticity  is  much 
higher  than  in  the  previous  two  cases.  Also,  the  structures  appear  to  be  much  thinner  and 
more  elongated  in  the  axial  direction.  Figure  5.65  shows  two  long  axial  vortices  which  are  much 
stronger  and  much  thinner  than  for  the  cases  with  Moo  =  0.2  and  Moo  =  1-2. 

Flow  visualization  using  particles  as  shown  in  figures  5.69  through  5.71  was  performed  similar 
to  the  subsonic  case.  As  before,  the  particles  were  introduced  at  a  ring  near  the  corner  of  the 
base  {x  =  0.025,  r  =  0.475).  The  plots  do  not  indicate  the  structures  as  clearly  as  in  the  cases 
for  Moo  =  0.2  and  Moo  =  1-2.  The  cross-sectional  cuts,  however,  indicate  that  the  longitudinal 
vortices  are  more  confined  in  the  shear  layer  than  for  the  cases  with  Moo  =  0.2  and  Moo  =  1-2. 
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Figure  5.53:  Isolines  of  instantaneous  total  vorticity  for  Reu  =  30,000  and  Moo  =  2.46  at 
t  =  178.1,  x-z  plane  (top)  and  x-y  plane  (bottom)  (filtered  DNS). 
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Figure  5.57:  Isosurfaces  of  instantaneous  total  vorticity  (|a;|  =  7)  for  Mq©  =  2.46,  I^d  =  30,000 
at  t  =  174.35  (filtered  DNS). 


Figure  5.58:  Isosurfaces  of  instantaneous  total  vorticity  (|n>|  =  7)  for  Moo  =  2.46,  ifep  =  30,000 
at  t  =  176.35  (filtered  DNS). 


Figure  5.59:  Isosurfaces  of  instantaneous  total  vorticity  (|c<;|  =  7)  for  M^o  =  2.46,  ifei)  =  30,000 
at  t  =  178.1  (filtered  DNS). 
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Figure  5.60:  Isosurfaces  of  instantaneous  azimuthal  vorticity  (|a;e|  =  7)  for  Moo  =  2.46,  = 

30,000  at  t  =  174.35  (filtered  DNS). 


Figure  5.61:  Isosurfaces  of  instantaneous  azimuthal  vorticity  (|uie|  =  7)  for  Moo  =  2.46,  ifej?  = 
30,000  at  t  =  176.35  (filtered  DNS). 


Figure  5.62:  Isosurfaces  of  instantaneous  azimuthal  vorticity  (loi^l  =  7)  for  Moo  —  2.46,  ife/j  = 
30,000  at  t  =  178.1  (filtered  DNS). 


Figure  5.63:  Isosurfaces  of  instantaneous  axial  vorticity  (|a;x|  =  4)  for  Moo  =  2.46,  Reo  =  30, 000 
at  t  =  174.35  (filtered  DNS). 
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Figure  5.64:  Isosurfaces  of  instantaneous  axial  vorticity  (Iwil  =  4)  for  Moo  =  2.46,  Reo  =  30, 000 
at  t  =  176.35  (filtered  DNS). 


Figure  5.65:  Isosurfaces  of  instantaneous  axial  vorticity  (|a;x|  =  4)  for  Mq©  =  2.46,  ife£)  =  30, 000 
at  t  =  178.1  (filtered  DNS). 


Figure  5.66:  Isosurfaces  of  instantaneous  pressure  deviation  firom  the  time- averaged  mean  flow 
(Ap  =  -0.02)  for  Moo  =  2.46,  I^d  =  30,000  at  t  =  174.35  (filtered  DNS). 


Figure  5.67:  Isosiurfaces  of  instantaneous  pressure  deviation  from  the  time-averaged  mean  flow 
(Ap  =  -0.02)  for  Moo  =  2.46,  Rsd  =  30,000  at  t  =  176.35  (filtered  DNS). 
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Figure  5.68:  Isosurfaces  of  instantaneous  pressure  deviation  from  the  time-averaged  mean  flow 
(Ap  =  —0.02)  for  Moo  =  2.46,  =  30,000  at  t  =  178.1  (Altered  DNS). 


Figure  5.69:  Flow  visualization  of  instantaneous  flow  field  via  particles  for  Moo  =  2.46,  ifeo  = 
30,000,  viewed  from  6  =  7r/2  (top)  and  from  0  =  0  (bottom)  (particles  introduced  at  a:  = 
0.025, r  =  0.475;  filtered  DNS). 

The  previous  cases  for  Moo  —  0.2  and  Moo  =  1-2  have  shown  that  the  unsteady  structures 
have  a  strong  influence  on  the  global  time-averaged  flow  field.  For  the  Moo  =  2.46  case,  the 
time-averaged  flow  is  shown  in  figures  5.72  and  5.73  in  the  form  of  isolines  of  the  pressure  and 
density.  As  before,  the  plots  exhibit  a  drastic  diff'erence  between  the  time-averaged  flow  and 
the  axisymmetric  flow  (compare  figures  4.8  and  4,9).  The  pressure  distribution  in  the  radial 
direction  for  the  time-averaged  flow  is  almost  constant  near  the  base.  This  is  in  agreement  with 
experimental  observations. 

Figure  5.74  shows  the  pressure  distribution  along  the  base  normalized  by  the  static  pressure 
in  the  free  stream.  The  plot  confirms  that  this  is  due  to  the  presence  of  the  large  structures.  The 
base  pressure  distribution  is  almost  constant  over  the  radius.  Averaged  over  the  radius  the  base 
pressure  is  significantly  lower  than  for  the  axisymmetric  steady  calculation.  This  indicates,  that 
for  this  higher  Mach  number  the  dynamic  structures  also  have  a  large  influence  on  the  global 
flow  field  which  results  in  a  lower  base  pressure  and,  consequently,  in  a  higher  base  drag. 

In  summary,  the  calculations  indicate  that  the  flow  for  the  higher  Mach  number,  Moo  ~  2.46, 
exhibits  an  absolute  instability  at  a  global  Reynolds  number  of  Bep  =  30,000.  As  for  the 
subsonic  case  (Moo  =  0-2)  and  the  slightly  supersonic  case  (Moo  =  1-2)  the  absolute  instability 
leads  to  large  structures.  In  contrast  to  the  previous  cases,  the  structures  tend  to  break  down 
to  smaller  scales.  This  is  an  indication  that  the  flow  might  transition  to  turbulence  if  the 
calcuktions  were  better  resolved.  Therefore,  further  investigation  of  this  flow  was  performed 
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Figure  5.70:  Flow  visualization  of  instantaneous  flow  field  via  particles  for  Mqq  =  2.46,  Rsd  = 
30,000,  sheets  of  thickness  0.01  at  0  =  0  (top)  and  6  =  7r/2  (bottom),  (particles  introduced  at 
X  =  0.025, r  =  0.475;  filtered  DNS). 


employing  LES  (see  following  section). 
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Figure  5.71:  Flow  visualization  of  instantaneous  flow  field  via  particles  for  Moo  =  2.46,  ifez?  = 
30,000,  cross-sectional  cuts  of  thickness  0.01  at  a;  =  1.0,  x  =  2.0,  x  =  3.0,  and  x  =  4.0  (particles 
introduced  at  a:  =  0.025,  r  =  0.475;  filtered  DNS). 
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Figure  5.74:  Base  pressure  distributions  for  Moo  =  2.46  and  BeD  =  30, 000. 
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5.2  Large-Eddy  Simulations 

As  shown  in  section  5.1.3,  the  DNS  for  Mqo  =  2.46  and  Ren  =  30,000  without  spatial  filter¬ 
ing  terminated  after  the  local  velocity  gradients  near  the  stagnation  points  became  too  large. 
Therefore,  to  extend  the  calculations,  a  spatial  filter  was  applied  in  the  axial  and  the  radial 
directions.  The  results  shown  in  section  5.1.3  indicate  that  using  the  spatial  filter  (in  addition 
to  a  very  dissipative  second-order  filter  at  every  100th  time  step)  the  small  scale  structures  could 
be  dissipated.  However,  the  larger  structures  remained  in  the  flow  field  (see  section  5.1.3). 

The  application  of  a  second-order  filter  in  order  to  dissipate  the  small  scale  structures  is 
physically  not  very  meaningful.  In  order  to  investigate  the  true  behavior  of  the  flow  field  at  this 
Reynolds  number,  either  a  higher  resolution  for  the  DNS  is  necessary  or  a  turbulence  model  is 
required  to  model  the  unresolved  scales.  Since  the  goal  is  to  calculate  flow  fields  at  even  higher 
Reynolds  numbers  (up  to  0(10®)),  a  much  higher  resolution  would  be  required  for  the  DNS. 
However,  the  simulation  with  the  spatial  and  temporal  resolution  of  case  T7  (see  table  A.3) 
was  already  very  CPU  time  intensive.  A  typical  calculation  for  this  flow  field  required  about 
600  MBytes  of  main  memory  and  in  the  order  of  500  CPU  hours  on  a  single  processor  of  a 
CRAY  C90.  Even  more  memory  and  CPU  time  would  be  required  for  calculations  with  higher 
resolution  (about  5  GBytes  for  twice  the  resolution  in  physical  space). 

As  stated  in  the  introduction,  with  Reynolds-Averaged  Navier-Stokes  (RANS)  it  is  not  pos¬ 
sible  to  capture  the  relevant  physics  of  this  flow.  Employing  Large-Eddy  Simulation  (LES), 
on  the  other  hand,  will  allow  investigation  of  the  dynamics  of  the  large  structmes  and  their 
effect  on  the  local  flow  behavior.  The  implementation  of  LES  is  discussed  in  section  2.4.  For 
the  LES  the  same  sixth-order  acciurate  compact  difference  filter  was  applied  in  the  axial  and 
radial  directions  (see  also  section  5.1.3).  However,  no  additional  second  order  central  difference 
filter  was  applied.  The  Smagorinsky  constant  for  the  sub-grid  scale  model  was  chosen  to  be 
Cs  =  0.065.  This  is  the  same  value  that  has  been  used  for  incompressible  turbulent  boundary 
layer  calculations  [[Bachman  (1996)]].  It  is  not  clear  if  this  value  is  also  useful  for  the  present 
case.  In  addition,  a  constant  value  is  not  realistic  for  a  flow  with  such  complexity.  For  this  case, 
a  dynamic  sub-grid  scale  model  would  be  more  appropriate. 

In  the  following,  preliminary  results  are  presented  for  a  firee  stream  Mach  number  of  Mqo  = 
2.46  and  two  different  global  Reynolds  numbers,  Reu  —  30,000  and  Red  =  100,000.  The 
results  are  presented  in  form  of  isolines  of  instantaneous  total  vorticity,  isolines  of  axial  velocity, 
radial  velocity,  pressure,  density,  temperature,  local  Mach  number,  and  azimuthal  vorticity  for 
the  time-averaged  turbulent  flow  field,  and  isolines  of  root-mean-square  values  (RMS)  of  axial, 
radial  and  tangential  turbulent  intensities  as  well  as  axial-radial  Reynolds  stress  <  u'v'  >  and 
turbulent  kinetic  energy  (Here  <>  denotes  the  RMS  value). 

As  a  first  case,  the  LES  was  applied  to  the  previously  calculated  flow  field  at  Red  =  30, 000. 
Figure  5.75  shows  isolines  of  instantaneous  total  vorticity  in  the  flow  field.  A  direct  comparison 
with  figures  5.51  through  5.53  for  the  filtered  DNS  for  M^o  =  2.46  and  Red  =  30, 000  shows  that 
the  flow  field  for  the  LES  looks  very  similar  to  that  of  the  spatially  filtered  DNS.  However,  in  the 
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Figure  5.76:  Isolines  of  axial  velocity  for  the  time-averaged  flow  field  at  Moo  =  2.46  and  Rep  — 
30,000  (LES). 

DNS,  structures  appeared  in  the  free  shear  layer,  emanating  from  the  corner.  These  structures 
seem  to  be  suppressed  in  the  LES.  However,  structmres  are  noticeable  within  the  recirculating 
region  which  are  convected  upstream  to  the  base. 

As  for  the  time-averaged  spatially  filtered  DNS,  the  recirculating  region  of  the  time-averaged 
LES  is  much  shorter  {Lr  «  2.8)  than  for  the  steady  laminar  calculation  (Lr  w  3.5).  This  can  be 
seen  in  figure  5.76,  where  isolines  of  the  axial  velocity  of  the  time-averaged  turbulent  flow  field 
are  shown.  Figures  5.77  through  5.82  show  isolines  of  time-averaged  radial  velocity,  pressure, 
density,  temperature,  local  Mach  number,  and  aizimuthal  vorticity,  respectively.  Comparison 
with  the  results  for  the  steady  flow  calculation  reveals,  that  the  changes  in  the  flow  field  (for 
example  the  pressure  drop  downstream  of  the  corner)  are  much  more  intense  for  the  turbulent 
flow.  In  addition,  figure  5.78  shows  that,  similar  to  the  time-average  for  the  unsteady  DNS 
calculations  for  Moo  =  2.46  and  for  Moo  =  1-2,  the  pressure  distribution  along  the  base  is 
almost  constant.  This  is  in  good  agreement  with  experimental  observations. 


Figure  5.80:  Isolines  of  temperature  for  the  time-averaged  flow  fleld  at  Moo  =  2.46  and  Ben  = 
30,000  (LES). 


Figure  5.81:  Isolines  of  local  Mach  number  for  the  time-averaged  flow  field  at  Moo  =  2.46  and 
ifeD  =  30,000  (LES). 


Figure  5.82:  Isolines  of  azimuthal  vorticity  for  the  time-averaged  flow  field  at  Moo  =  2.46  and 
ifec  =  30,000  (LES). 
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For  further  comparison  with  experimental  results  from  Dutton  and  co-workers  [see,  for  ex¬ 
ample,  [Herrin  &:  Dutton  (1994c)]]  the  RMS  values  of  the  tmbulent  intensities,  the  Reynolds 
stress,  and  the  turbulent  kinetic  energy  are  shown  in  figures  5.83  through  5.87,  with 


axial  turbulence  intensity, 
radial  turbulence  intensity, 
tangential  turbulence  intensity, 
axial-radial  Reynolds  stress, 
turbulent  kinetic  energy. 


=  \/<  KK  >. 

ar  -  y/<  u'u'  >, 
ae  =  v'gv'g  >, 

^xr  =<  KK  >) 


The  values  of  the  turbulence  intensities,  the  Reynolds  stress  and  the  TKE  from  the  simulation 
are  slightly  below  the  experimental  measurements.  The  highest  level  of  turbulence  intensity 
in  the  simulations  can  be  found  in  the  free  shear  layer  just  upstream  of  the  stagnation  point 
(x  w  2.8).  This  finding  is  also  in  agreement  with  experimental  observations.  However,  the  results 
in  the  experimental  investigations  by  Dutton  and  co-workers  were  obtained  for  a  global  Reynolds 
number  Bsd  =  1, 600, 000  and  while  the  Reynolds  number  of  the  simulation  was  =  30, 000. 


Figure  5.83:  Isolines  of  axial  turbulence  intensity  at  Moo  =  2.46  and  Hsd  =  30,000  (LES). 


Figure  5.84:  Isolines  of  radial  turbulence  intensity  at  Mq©  =  2.46  and  Bsd  =  30,000  (LES). 
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Therefore,  other  simulations  were  carried  out,  where  the  global  Reynolds  number  was  raised 
to  =  100, 000  in  order  to  determine  the  effect  of  the  global  Reynolds  number  on  the  flow 
field  behavior.  Figure  5.88  shows  isolines  of  instantaneous  total  vorticity  for  Rep  =  100,000.  A 
comparison  with  figure  5.75  for  Rec  =  30,000  shows  that  more  small  scale  structures  appear 
to  be  present  for  the  higher  Reynolds  number  case  than  for  the  lower  Reynolds  number  case. 
The  time-averaged  flow  field,  which  is  shown  in  figures  5.89  through  5.95,  is  very  similar  to 
the  flow  field  at  the  lower  Reynolds  number,  except  that  the  gradients  in  the  flow  field  (for 
example  the  pressure  gradient  in  the  expansion  and  the  recompression  region)  are  more  intense 
ioT  Reu  =  100,000  . 

However,  at  the  higher  Reynolds  number  the  levels  of  the  RMS  values  of  the  turbulence 
intensities  are  almost  identical  with  the  levels  observed  in  experiments,  as  shown  in  figures  5.96 
through  5.100,  although  the  Reynolds  number  in  the  experiments  was  considerably  larger.  The 
averaged  flow  field  shows  a  somewhat  longer  recirculating  region  than  the  experimental  flow 
field.  This  could  be  a  consequence  of  stronger  instabilities  in  the  shear  layer  for  the  higher 
Reynolds  number  in  the  experiments.  On  the  other  hand,  the  discrepancy  could  also  be  due  to 
the  inadequate  subgrid-scale  model  (Smagorinsky,  with  constant  coefficient)  and/or  insufficient 
resolution  for  the  calculations  of  the  resolved  scale. 


Figure  5.88:  Isolines  of  instantaneous  total  vorticity  for  M^o  =  2.46  and  Rbd  =  100,000  (LES). 
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Figure  5.95:  Isolines  of  vorticity  for  the  time-averaged  flow  field  at  Moo  =  2.46  and  Ren  = 
100,000  (LES). 


Figure  5.96:  Isolines  of  axial  turbulence  intensity  at  Moo  =  2.46  and  Rep  =  100,000  (LES). 


Figure  5.98:  Isolines  of  tangential  turbulence  intensity  at  Moo  =  2.46  and  Rsp  =  100, 000  (LES). 
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Figure  5.100:  Isolines  of  turbulent  kinetic  energy  at  Moo  =  2.46  and  Reo  =  100,000  (LES). 
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Conclusions 


A  numerical  method  has  been  developed  for  studying  the  evolution  and  propagation  of  three- 
dimensional  structures  in  the  laminar  and  the  turbulent  wake  of  axisymmetric  bluff  bodies  with  a 
blunt  base  in  supersonic  flows.  The  method  is  based  on  the  complete  compressible  Navier-Stokes 
equations.  The  equations  are  solved  in  a  cylindrical  coordinate  system  using  finite  diflference 
approximations  of  fomth-order  accuracy  in  the  axial  and  radial  directions.  In  the  azimuthal 
direction  a  pseudo-spectral  method  is  employed.  The  method  is  explicit  in  time,  using  a  fourth- 
order  accurate  Runge-Kutta  scheme  for  the  time  integration. 

Axisymmetric  flows  have  been  calculated  by  solving  the  axisymmetric 
Navier-Stokes  equations.  All  axisymmetric  flow  fields  have  been  found  to  converge  to  a  steady 
state.  These  flow  fields  are  then  used  as  initial  conditions  for  the  time-dependent  three- 
dimensional  calculations.  The  evolution  of  three-dimensional  disturbances  was  investigated  by 
distmrbing  the  axisymmetric  flow  either  through  a  blowing  and  suction  slot  near  the  base  of  the 
body,  or  by  introducing  local  single  pulse  disturbances  in  the  flow  field. 

For  code  validation,  comparison  of  the  results  for  subsonic  calculations  was  made  with  wa¬ 
ter  channel  experiments  and  incompressible  simulations.  It  was  found  that  for  the  absolutely 
unstable  flow  at  a  Reynolds  number  ifep  =  1,000  the  compressible  calculations  show  similar 
structures  as  observed  in  the  experiments  and  the  incompressible  calculations.  For  a  continu¬ 
ously  excited  flow  field,  amplitude  distributions  of  the  disturbances  in  the  wake  region  closely 
matched  those  of  the  water  channel  experiments. 

Direct  Numerical  Simulations  (DNS)  were  carried  out  for  three  different  free-stream  Mach 
numbers.  Absolute  instabilities  were  observed  for  the  following  Mach  number  -  Reynolds  number 
combinations:  Mqo  =  0.2  and  Bej)  —  1,000,  Moo  =  1-2  and  ifep  =  4,000,  and  Moo  =  2.46  and 
Bed  =  30,000.  For  Moo  —  2.46,  no  absolute  instability  was  found  at  lower  Reynolds  numbers 
{Bed  <  30,000).  Rather,  the  flow  always  returned  to  the  initial  axisymmetric  steady  state.  For 
the  three  Mach  number  and  Reynolds  number  combinations  shown  above,  however,  the  flow 
deviated  from  the  axisymmetric  steady  state  and  reached  a  state  with  non-periodic  fluctuations. 
At  this  point,  structures  appeared  in  the  flow  field,  originating  near  the  corner  of  the  blunt  base 
and  increasing  in  strength  as  they  were  convected  downstream.  For  all  Mach  number  -  Reynolds 
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number  combinations  the  structures  consist  of  a  region  of  high  azimuthal  vorticity  followed  by 
two  regions  of  axial  vorticity  of  opposite  signs  (indicating  two  counter-rotating  vortices).  For  the 
cases  Moo  =  0.2  and  Moo  =  1-2  the  observed  structures  can  be  explained  by  the  superposition 
of  two  helical  modes  (A:  =  1  and  k  =  —1).  For  the  higher  Mach  number  case  {Moo  =  2.46), 
however,  the  vortices  are  more  confined  within  the  shear  layer. 

The  time  dependent  large  structures  have  a  strong  influence  on  the  global,  time-averaged 
flow.  This  can  been  shown  by  comparing  the  time  average  of  the  unsteady  three-dimensional 
flow  with  the  initial  axisymmetric  steady  state.  For  the  time-averaged  flow  the  base  pressure 
distribution  is  almost  constant  over  the  radius.  In  addition,  the  integral  of  the  base  pressure  over 
the  entire  base  area  is  much  lower  for  the  time-averaged  flow  than  for  the  steady  axisymmetric 
calculations.  This  indicates  that  an  additional  contribution  of  base  drag  is  added  due  to  the 
dynamical  behavior  of  the  large  structmes. 

For  the  case  of  Moo  =  2.46  and  Rej)  =  30,000  attempts  of  DNS  were  terminated  because  of 
lack  of  resolution  in  the  calculation.  Therefore,  Large-Eddy  Simulations  (LES)  were  employed  for 
this  and  higher  Reynolds  numbers.  The  results  of  the  LES  show  that  dominant  large  structures 
exist  in  the  flow,  rendering  the  flow  fleld  highly  unsteady.  However,  the  evolution  of  large 
oscillating  structures  in  the  shear  layer  emanating  from  the  corner  appears  to  be  suppressed, 
possibly  due  to  the  inadequacies  of  the  subgrid-scale  model  used  in  the  LES  and/or  due  to 
insufiicient  resolution  for  the  calculation  of  the  resolved  scales.  For  LES  at  a  higher  global 
Reynolds  number  {Reo  =  100, 000)  turbulence  levels  are  reached  that  closely  match  those  of  the 
experiments  for  much  larger  Reynolds  numbers  {Ren  =  1,600,000)  [see,  for  example,  Herrin  & 
Dutton  (1994c)]. 

However,  other  time-averaged  quantities  do  not  agree  well  with  the  experimental  flndings. 
The  length  of  the  recirculating  region  for  the  LES  calculations  is  Lr  «  2.8,  whereas  in  the 
experimental  investigation  it  was  found  to  be  Lr  w  1.3  [Herrin  &  Dutton  (1994c)].  In  addition, 
the  levels  of  the  kinetic  energy  and  of  the  Reynolds  shear  stress  in  the  experiments  reveal  that 
the  shear  layer  is  highly  unstable  immediately  downstream  of  the  corner.  As  mentioned  above, 
this  instability  was  also  observed  in  the  preliminary  DNS  for  Rsd  =  30,000,  but  not  in  the 
LES  calculations.  Therefore,  further  investigations  are  needed  to  explore  the  influence  of  the 
subgrid-scale  model  and/or  the  resolution  used  in  the  calculations. 

In  general,  the  results  of  the  simulations  indicate  that  it  may  be  possible  to  influence  the 
base  pressure  (i.e.  the  total  aerodynamic  drag)  of  axisymmetric  bodies  with  a  blunt  base  by 
influencing  the  generation  and/or  the  dynamics  of  the  large  structures.  This  can  be  done  either 
by  passive  means  (base  bleed,  boat-tailing,  etc.)  which  is  cmrently  investigated  by  Dutton  and 
co-workers,  or  by  active  means  (introducing  controlled  periodic  disturbances). 
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TABLES  OF  PARAMETERS 


Case 

Schwarz 

A1 

A2 

A3 

A4 

A5 

A6 

Moo 

N/A 

0.2 

2.46 

1.2 

Bed 

1,000 

30,000 

4,000 

Poo[bar] 

N/A 

1.01325 

0.3207851 

Toc[K] 

N/A 

72.0 

133.0 

N/A 

1,004.9157 

Pr 

N/A 

0.7 

7 

N/A 

1.4 

N 

209 

112 

223 

112 

640 

1279 

420 

M 

129 

50 

99 

50 

120 

239 

400 

Xl 

-1.5 

-0.71956 

-1.5625 

Xn 

5.0 

18.15132 

14.48345 

5.0 

tm 

4.0 

6.2841 

6.0435 

4.0 

Ae 

0.03125 

0.025 

0.0125 

0.025 

0.004427 

0.0022135 

0.015625 

At] 

0.03125 

0.025 

0.0125 

0.025 

0.01 

0.005 

0.01 

Sx 

1.0 

0.4391 

0.7 

1.0 

1.0 

0.4 

0.3 

1.0 

At 

N/A 

0.0014 

0.00068 

0.0014 

0.0019 

0.00096 

0.0019 

Po 

N/A 

1.0 

N/A 

■Aq 

N/A 

0.001 

N/A 

Table  A.l:  Parameters  for  axisymmetric  calculations. 
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T1 

T2 

T3 

T4 

T5 

Moo 

0.2 

Ren 

1,000 

Poa[bar] 

1.01325 

Too[K] 

72.0 

^  kqK 

1004.9157 

Pr 

0.7 

7 

1.4 

N 

112 

223 

112 

M 

50 

99 

50 

Xl 

-1.5 

XN 

18.15132 

rM 

6.2841 

0.025 

0.0125 

0.025 

At) 

0.025 

0.0125 

0.025 

Sx 

0.4391 

Sr 

0.4 

At 

0.0014 

0.00068 

0.0014 

fcmax 

4 

8 

S^D 

0.0318 

At-Pulse 

0.159 

Ampl. 

0.001 

0.1 

Table  A.2:  Parameters  for  three-dimensional  subsonic  DNS  calculations. 
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T6 

T7 

Moo 

1.2 

2.46 

Ren 

4,000 

30,000 

Pcx,[bar] 

0.3207851 

Too[K] 

133.0 

^  kgK 

1004.9157 

Pr 

0.7 

7 

1.4 

N 

420 

640 

M 

400 

120 

Xl 

-0.71956 

-1.5625 

XN 

5.0 

14.48345 

tm 

4.0 

6.0435 

0.015625 

0.004427 

At] 

0.01 

0.01 

Sx 

1.0 

0.7 

Sf 

1.0 

0.3 

At 

0.0019 

0.0019 

fcmax 

4 

8 

StD 

At-Pulse 

Ampl. 

0.0001 

Table  A.3:  Parameters  for  three-dimensional  supersonic  DNS  calculations. 
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FREE  STREAM  BOUNDARY 
CONDITIONS 


The  boundary  conditions  at  the  free-stream  boundary  are  discussed  in  detail  in  Harris(1995). 
The  equations  change  only  slightly  for  the  axisymmetric  coordinate  system  used  here.  Depend- 


ing  on  the  free  stream  Mach  number 

(Moo),  three  different  kinds  of  boundary  conditions  are 

considered: 

Moo  <  1  : 

subsonic  boundary  conditions. 

(B.l) 

Moo  =  1-2  : 

characteristic  boundary  conditions. 

(B.2) 

Moo  =  2.46  : 

Thompson  boundary  conditions. 

(B.3) 

The  three  different  methods  are  described  in  detail  in  the  next  sections. 


B.l  Subsonic  Boundary  Conditions 

For  subsonic  flow  fields,  where  the  local  Mach  number  at  the  free  stream  is  always  less  than 
one,  the  flow  field  is  assumed  to  have  reached  nearly  axisymmetric  state.  Therefore,  derivatives 
in  radial  direction  of  any  flow  variable  are  neglected.  This  leads  to  the  following  boundary 
conditions: 

|:  W  =  0.  (B.4) 

Prom  continuity  it  follows  that 

^  [rvr]  =  0.  (B.5) 

B.2  Characteristic  Boundary  Conditions 

The  second  free-stream  boundary  condition  (referred  to  hereafter  as  the  Characteristic  Boundary 
Condition)  is  applied  when  the  free  stream  Mach  number  is  Mqo  =  1-2.  For  steady,  isentropic 
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two-dimensional  flow, 

1  =  0  (B.6) 

is  used  for  all  variables,  where  c  is  the  outward-travelling  characteristic  direction  along  the 
Mach  angle.  Corrected  for  flow  direction,  the  angle  of  the  characteristic  is  given  at  the  upper 
free  stream  by 

where  the  local  Mach  number  is 


The  lo-terms  are  not  included  because  the  flow  is  presumably  (for  each  /-plane)  steady  and 
axisymmetric,  thus  w  =  0.  The  values  of  u,v,w,T,  and  p  are  thus  found  using 

(B  9) 

where  {idiJ  —  1)  and  (ic2,J  —  2)  are  backward  locations  along  c,  starting  from  (/,  J).  The 
values  at  those  locations  are  found  by  linear  interpolation  of  the  values  at  surrounding  points, 
e.g.,  with  ica  <  id  <  icb, 

=  (1  _  0:)^*“’'^“^  +  (B.IO) 

B.3  Thompson  Boundary  Conditions 

The  third  free-stream  boundary  conditions  used  were  suggested  by  Thompson  [1987].  Here,  the 
integration  is  extended  to  the  boundary  itself,  where  the  integration  in  the  radial  direction  is 
done  in  a  manner  which  should  allow  distiurbances  to  move  out  of,  but  not  into  the  domain. 

The  method  is  described  in  detail  in  Harris  (1995).  It  is  based  on  the  fact  that  information 
can  travel  along  the  outgoing  characteristic,  but  is  not  allowed  to  propagate  along  the  incoming 
characteristic.  Thus,  only  outgoing  information  is  allowed  at  the  free  stream  boundary.  There¬ 
fore,  a  special  treatment  is  required  for  the  convective  terms  in  the  radial  direction,  while  for 
the  viscous  terms  one-sided  difi'erences  are  applied  for  derivatives  with  respect  to  r. 

This  condition  was  found  to  lead  to  reflections  of  waves  for  subsonic  and  slightly  supersonic 
flows.  Only  for  a  free  stream  Mach  number  of  M^o  =  2.46  no  reflections  could  be  observed. 
For  this  reason,  this  boundary  condition  was  only  applied  for  that  particular  free  stream  Mach 
number.  For  the  Moo  =  1-2  case  the  Characteristic  Boundary  Conditions  were  employed  (see 
section  B.2). 


Appendix  C 


BOUNDARY  CONDITIONS  AT 
THE  AXIS  OF  SYMMETRY 


The  boundary  conditions  at  the  symmetry  axis  are  non-physical  and  appear  in  the  numerical 
scheme  because  of  the  Fourier  decomposition  of  the  flow  variables  in  circumferential  direction 
{6).  The  conditions  are  independent  of  time.  Therefore  the  time  indices  are  dropped  in  the 
following  derivation. 

The  derivation  of  the  boundary  condition  is  similar  to  the  one  for  incompressible  computation 
as  shown  for  example  in  Tourbier  (1991).  For  the  cylindrical  coordinate  system  (a;,  r,  9)  the  axis 
(r  =  0)  forms  a  singularity,  that  needs  special  treatment  in  the  numerical  scheme.  This  results 
mainly  in  the  treatment  of  vectors,  that  are  non-parallel  to  the  axis  and  terms  with  division  by 
the  radial  coordinate  (r). 
a)  Velocity  Vector  at  the  Axis 


Figure  C.l:  Decomposition  of  the  velocity  vector  within  the  field. 

The  velocity  vector  (vx,  Vr,ve)  is  related  to  the  cartesian  velocity  vector  (ux,  Vy,  v^)  as  follows: 
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or 


Vx  —  Ui, 

Vr  —  VyCosO  +  VzSinO,  (C.2) 

ve  =  —Vy  sin  B  +  Vz  cos  6,  (C.3) 

Vx  =  Ux,  (C.4) 

Vy  =  Vr  cos  B  —  Vff  sin  0,  (C.5) 

Vz  =  Vr  sin  9  +  ve  cos  9.  (C.6) 


The  Fourier  decomposition  of  the  velocity  vector  components  yields 

K 

Vx(x,r,9)  =  (C.7) 

k=-K 

K 

Vr{x,r,9)  =  53  '»r,k{x,r)e^’^\  (C.8) 

k~-K 

K 

ve{x,T,9)  =  53  .  (C.9) 

k=-K 


Figure  C.2:  Decomposition  of  the  velocity  vector  at  the  axis  of  symmetry. 

At  the  axis  (r  =  0)  all  angles  (0)  belong  physically  to  the  same  point  {y,  z)  and  to  unique 
values  of  Vy  and  Vz.  In  addition,  all  scalar  values  as  well  as  the  axial  velocity  component  (vx) 
are  the  same  for  all  angles  (9).  Therefore  for  all  scalar  fields  and  the  axial  velocity  component: 


po(a;,0)  ^  0, 
To(a:,0)  ^  0, 


(C.IO) 

(C.ll) 


and 


nONS  AT  THE  AXIS  OF  SYMMETRY 

13'4 

Po(a;,0)  7^0, 

(C.12) 

eo(a:,0)  7^0, 

(C.13) 

^  0) 

(C.14) 

Pfc(a;,0)  =  0, 

(C.15) 

fkix,0)  =  0, 

(C.16) 

Pik(«,0)  =  0, 

(C.17) 

efe(x,0)  =  0, 

(C.18) 

vi,fe(«,0)  =  0, 

(C.19) 

for  A:  ^  0. 

The  radial  and  the  azimuthal  component  of  the  velocity  vector  at  the  axis  are  represented 
by 


K 

Ur  (®)  0, 0)  =  ^  Vr,kix,  0)  [cos(fc0)  + 1  sin(A:0)] ,  (C.20) 

k--K 
K 

Vg{x,0,6)  =  Ufl,*,(a:,0)  [cos(A:0) +isin(A:0)] .  (C.21) 

k=-K 

Since  the  functions  sin(A:0)  and  cos(A:0)  build  an  orthogonal  basis,  it  follows,  that  only  the 
Fourier  components  k  —  1  and  A:  =  —  1  can  be  non-zero  when  comparing  the  components  to 
equation  C.5.  Therefore, 


Ur(a;,O,0)  =  [up,i(a;,0) -l-Ur,-i(®»O)]cos0[up,i(a:,O) -i;rf-i(x,O)]sin0,  (C.22) 

vg{x,0,6)  =  [u0,i(a;,O)  H-V0"_i(a:,O)]cos0[u0,i(a:,O) -fU0"_i(a;,O)]cos0.  (C.23) 

Combining  these  results  leads  to  the  final  conditions  for  the  radial  and  the  azimuthal  velocity 
components  at  the  symmetry  axis: 


v;,i(a:,  0)  +  ve,i{x,  0)  =  0,  (C.24) 

Vr'-i  (a:,  0)  -  ve~-i  {x,  0)  =  0.  (C.25) 

b)  Treatment  of  the  Division  by  r  Near  the  Axis 

For  the  division  by  r  at  the  axis  of  symmetry  the  rule  of  I’Hopital  was  applied.  For  the 
example  of  the  term  ^  this  leads  to 


lim  — = 
r-»0  r 


dVr 

dr 


This  rule  can  be  applied,  as  long  as  the  numerator  vanishes  at  r  =  0. 


(C.26) 
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